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.
25 const BASE alpha_real = CONST_REAL0(alpha);
26 const BASE alpha_imag = CONST_IMAG0(alpha);
28 const BASE beta_real = CONST_REAL0(beta);
29 const BASE beta_imag = CONST_IMAG0(beta);
31 if ((alpha_real == 0.0 && alpha_imag == 0.0)
32 && (beta_real == 1.0 && beta_imag == 0.0))
35 if (Order == CblasRowMajor) {
43 uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;
44 side = (Side == CblasLeft) ? CblasRight : CblasLeft;
47 /* form y := beta*y */
48 if (beta_real == 0.0 && beta_imag == 0.0) {
49 for (i = 0; i < n1; i++) {
50 for (j = 0; j < n2; j++) {
51 REAL(C, ldc * i + j) = 0.0;
52 IMAG(C, ldc * i + j) = 0.0;
55 } else if (!(beta_real == 1.0 && beta_imag == 0.0)) {
56 for (i = 0; i < n1; i++) {
57 for (j = 0; j < n2; j++) {
58 const BASE Cij_real = REAL(C, ldc * i + j);
59 const BASE Cij_imag = IMAG(C, ldc * i + j);
60 REAL(C, ldc * i + j) = beta_real * Cij_real - beta_imag * Cij_imag;
61 IMAG(C, ldc * i + j) = beta_real * Cij_imag + beta_imag * Cij_real;
66 if (alpha_real == 0.0 && alpha_imag == 0.0)
69 if (side == CblasLeft && uplo == CblasUpper) {
71 /* form C := alpha*A*B + C */
73 for (i = 0; i < n1; i++) {
74 for (j = 0; j < n2; j++) {
75 const BASE Bij_real = CONST_REAL(B, ldb * i + j);
76 const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
77 const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
78 const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
79 BASE temp2_real = 0.0;
80 BASE temp2_imag = 0.0;
82 const BASE Aii_real = CONST_REAL(A, i * lda + i);
83 /* const BASE Aii_imag = 0.0; */
84 REAL(C, i * ldc + j) += temp1_real * Aii_real;
85 IMAG(C, i * ldc + j) += temp1_imag * Aii_real;
87 for (k = i + 1; k < n1; k++) {
88 const BASE Aik_real = CONST_REAL(A, i * lda + k);
89 const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
90 const BASE Bkj_real = CONST_REAL(B, ldb * k + j);
91 const BASE Bkj_imag = CONST_IMAG(B, ldb * k + j);
92 REAL(C, k * ldc + j) += Aik_real * temp1_real - (-Aik_imag) * temp1_imag;
93 IMAG(C, k * ldc + j) += Aik_real * temp1_imag + (-Aik_imag) * temp1_real;
94 temp2_real += Aik_real * Bkj_real - Aik_imag * Bkj_imag;
95 temp2_imag += Aik_real * Bkj_imag + Aik_imag * Bkj_real;
97 REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
98 IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
102 } else if (side == CblasLeft && uplo == CblasLower) {
104 /* form C := alpha*A*B + C */
106 for (i = 0; i < n1; i++) {
107 for (j = 0; j < n2; j++) {
108 const BASE Bij_real = CONST_REAL(B, ldb * i + j);
109 const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
110 const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
111 const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
112 BASE temp2_real = 0.0;
113 BASE temp2_imag = 0.0;
114 for (k = 0; k < i; k++) {
115 const BASE Aik_real = CONST_REAL(A, i * lda + k);
116 const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
117 const BASE Bkj_real = CONST_REAL(B, ldb * k + j);
118 const BASE Bkj_imag = CONST_IMAG(B, ldb * k + j);
119 REAL(C, k * ldc + j) += Aik_real * temp1_real - (-Aik_imag) * temp1_imag;
120 IMAG(C, k * ldc + j) += Aik_real * temp1_imag + (-Aik_imag) * temp1_real;
121 temp2_real += Aik_real * Bkj_real - Aik_imag * Bkj_imag;
122 temp2_imag += Aik_real * Bkj_imag + Aik_imag * Bkj_real;
125 const BASE Aii_real = CONST_REAL(A, i * lda + i);
126 /* const BASE Aii_imag = 0.0; */
127 REAL(C, i * ldc + j) += temp1_real * Aii_real;
128 IMAG(C, i * ldc + j) += temp1_imag * Aii_real;
130 REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
131 IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
135 } else if (side == CblasRight && uplo == CblasUpper) {
137 /* form C := alpha*B*A + C */
139 for (i = 0; i < n1; i++) {
140 for (j = 0; j < n2; j++) {
141 const BASE Bij_real = CONST_REAL(B, ldb * i + j);
142 const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
143 const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
144 const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
145 BASE temp2_real = 0.0;
146 BASE temp2_imag = 0.0;
148 const BASE Ajj_real = CONST_REAL(A, j * lda + j);
149 /* const BASE Ajj_imag = 0.0; */
150 REAL(C, i * ldc + j) += temp1_real * Ajj_real;
151 IMAG(C, i * ldc + j) += temp1_imag * Ajj_real;
153 for (k = j + 1; k < n2; k++) {
154 const BASE Ajk_real = CONST_REAL(A, j * lda + k);
155 const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
156 const BASE Bik_real = CONST_REAL(B, ldb * i + k);
157 const BASE Bik_imag = CONST_IMAG(B, ldb * i + k);
158 REAL(C, i * ldc + k) += temp1_real * Ajk_real - temp1_imag * Ajk_imag;
159 IMAG(C, i * ldc + k) += temp1_real * Ajk_imag + temp1_imag * Ajk_real;
160 temp2_real += Bik_real * Ajk_real - Bik_imag * (-Ajk_imag);
161 temp2_imag += Bik_real * (-Ajk_imag) + Bik_imag * Ajk_real;
163 REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
164 IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
168 } else if (side == CblasRight && uplo == CblasLower) {
170 /* form C := alpha*B*A + C */
172 for (i = 0; i < n1; i++) {
173 for (j = 0; j < n2; j++) {
174 const BASE Bij_real = CONST_REAL(B, ldb * i + j);
175 const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
176 const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
177 const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
178 BASE temp2_real = 0.0;
179 BASE temp2_imag = 0.0;
180 for (k = 0; k < j; k++) {
181 const BASE Ajk_real = CONST_REAL(A, j * lda + k);
182 const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
183 const BASE Bik_real = CONST_REAL(B, ldb * i + k);
184 const BASE Bik_imag = CONST_IMAG(B, ldb * i + k);
185 REAL(C, i * ldc + k) += temp1_real * Ajk_real - temp1_imag * Ajk_imag;
186 IMAG(C, i * ldc + k) += temp1_real * Ajk_imag + temp1_imag * Ajk_real;
187 temp2_real += Bik_real * Ajk_real - Bik_imag * (-Ajk_imag);
188 temp2_imag += Bik_real * (-Ajk_imag) + Bik_imag * Ajk_real;
191 const BASE Ajj_real = CONST_REAL(A, j * lda + j);
192 /* const BASE Ajj_imag = 0.0; */
193 REAL(C, i * ldc + j) += temp1_real * Ajj_real;
194 IMAG(C, i * ldc + j) += temp1_imag * Ajj_real;
196 REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
197 IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
202 BLAS_ERROR("unrecognized operation");