1 /* blas/source_syrk_r.h
3 * Copyright (C) 2001, 2007 Brian Gough
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 const BASE alpha_real = CONST_REAL0(alpha);
25 const BASE alpha_imag = CONST_IMAG0(alpha);
27 const BASE beta_real = CONST_REAL0(beta);
28 const BASE beta_imag = CONST_IMAG0(beta);
30 if ((alpha_real == 0.0 && alpha_imag == 0.0)
31 && (beta_real == 1.0 && beta_imag == 0.0))
34 if (Order == CblasRowMajor) {
36 /* FIXME: original blas does not make distinction between Trans and ConjTrans?? */
37 trans = (Trans == CblasNoTrans) ? CblasNoTrans : CblasTrans;
39 uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;
40 trans = (Trans == CblasNoTrans) ? CblasTrans : CblasNoTrans;
43 /* form y := beta*y */
44 if (beta_real == 0.0 && beta_imag == 0.0) {
45 if (uplo == CblasUpper) {
46 for (i = 0; i < N; i++) {
47 for (j = i; j < N; j++) {
48 REAL(C, ldc * i + j) = 0.0;
49 IMAG(C, ldc * i + j) = 0.0;
53 for (i = 0; i < N; i++) {
54 for (j = 0; j <= i; j++) {
55 REAL(C, ldc * i + j) = 0.0;
56 IMAG(C, ldc * i + j) = 0.0;
60 } else if (!(beta_real == 1.0 && beta_imag == 0.0)) {
61 if (uplo == CblasUpper) {
62 for (i = 0; i < N; i++) {
63 for (j = i; j < N; j++) {
64 const BASE Cij_real = REAL(C, ldc * i + j);
65 const BASE Cij_imag = IMAG(C, ldc * i + j);
66 REAL(C, ldc * i + j) = beta_real * Cij_real - beta_imag * Cij_imag;
67 IMAG(C, ldc * i + j) = beta_real * Cij_imag + beta_imag * Cij_real;
71 for (i = 0; i < N; i++) {
72 for (j = 0; j <= i; j++) {
73 const BASE Cij_real = REAL(C, ldc * i + j);
74 const BASE Cij_imag = IMAG(C, ldc * i + j);
75 REAL(C, ldc * i + j) = beta_real * Cij_real - beta_imag * Cij_imag;
76 IMAG(C, ldc * i + j) = beta_real * Cij_imag + beta_imag * Cij_real;
82 if (alpha_real == 0.0 && alpha_imag == 0.0)
85 if (uplo == CblasUpper && trans == CblasNoTrans) {
87 for (i = 0; i < N; i++) {
88 for (j = i; j < N; j++) {
91 for (k = 0; k < K; k++) {
92 const BASE Aik_real = CONST_REAL(A, i * lda + k);
93 const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
94 const BASE Ajk_real = CONST_REAL(A, j * lda + k);
95 const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
96 temp_real += Aik_real * Ajk_real - Aik_imag * Ajk_imag;
97 temp_imag += Aik_real * Ajk_imag + Aik_imag * Ajk_real;
99 REAL(C, i * ldc + j) += alpha_real * temp_real - alpha_imag * temp_imag;
100 IMAG(C, i * ldc + j) += alpha_real * temp_imag + alpha_imag * temp_real;
104 } else if (uplo == CblasUpper && trans == CblasTrans) {
106 for (i = 0; i < N; i++) {
107 for (j = i; j < N; j++) {
108 BASE temp_real = 0.0;
109 BASE temp_imag = 0.0;
110 for (k = 0; k < K; k++) {
111 const BASE Aki_real = CONST_REAL(A, k * lda + i);
112 const BASE Aki_imag = CONST_IMAG(A, k * lda + i);
113 const BASE Akj_real = CONST_REAL(A, k * lda + j);
114 const BASE Akj_imag = CONST_IMAG(A, k * lda + j);
115 temp_real += Aki_real * Akj_real - Aki_imag * Akj_imag;
116 temp_imag += Aki_real * Akj_imag + Aki_imag * Akj_real;
118 REAL(C, i * ldc + j) += alpha_real * temp_real - alpha_imag * temp_imag;
119 IMAG(C, i * ldc + j) += alpha_real * temp_imag + alpha_imag * temp_real;
123 } else if (uplo == CblasLower && trans == CblasNoTrans) {
125 for (i = 0; i < N; i++) {
126 for (j = 0; j <= i; j++) {
127 BASE temp_real = 0.0;
128 BASE temp_imag = 0.0;
129 for (k = 0; k < K; k++) {
130 const BASE Aik_real = CONST_REAL(A, i * lda + k);
131 const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
132 const BASE Ajk_real = CONST_REAL(A, j * lda + k);
133 const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
134 temp_real += Aik_real * Ajk_real - Aik_imag * Ajk_imag;
135 temp_imag += Aik_real * Ajk_imag + Aik_imag * Ajk_real;
137 REAL(C, i * ldc + j) += alpha_real * temp_real - alpha_imag * temp_imag;
138 IMAG(C, i * ldc + j) += alpha_real * temp_imag + alpha_imag * temp_real;
142 } else if (uplo == CblasLower && trans == CblasTrans) {
144 for (i = 0; i < N; i++) {
145 for (j = 0; j <= i; j++) {
146 BASE temp_real = 0.0;
147 BASE temp_imag = 0.0;
148 for (k = 0; k < K; k++) {
149 const BASE Aki_real = CONST_REAL(A, k * lda + i);
150 const BASE Aki_imag = CONST_IMAG(A, k * lda + i);
151 const BASE Akj_real = CONST_REAL(A, k * lda + j);
152 const BASE Akj_imag = CONST_IMAG(A, k * lda + j);
153 temp_real += Aki_real * Akj_real - Aki_imag * Akj_imag;
154 temp_imag += Aki_real * Akj_imag + Aki_imag * Akj_real;
156 REAL(C, i * ldc + j) += alpha_real * temp_real - alpha_imag * temp_imag;
157 IMAG(C, i * ldc + j) += alpha_real * temp_imag + alpha_imag * temp_real;
162 BLAS_ERROR("unrecognized operation");