1 /* blas/source_trmm_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.
23 const int nonunit = (Diag == CblasNonUnit);
24 int side, uplo, trans;
26 if (Order == CblasRowMajor) {
31 trans = (TransA == CblasConjTrans) ? CblasTrans : TransA;
35 side = (Side == CblasLeft) ? CblasRight : CblasLeft;
36 uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;
37 trans = (TransA == CblasConjTrans) ? CblasTrans : TransA;
40 if (side == CblasLeft && uplo == CblasUpper && trans == CblasNoTrans) {
42 /* form B := alpha * TriU(A)*B */
44 for (i = 0; i < n1; i++) {
45 for (j = 0; j < n2; j++) {
49 temp = A[i * lda + i] * B[i * ldb + j];
51 temp = B[i * ldb + j];
54 for (k = i + 1; k < n1; k++) {
55 temp += A[lda * i + k] * B[k * ldb + j];
58 B[ldb * i + j] = alpha * temp;
62 } else if (side == CblasLeft && uplo == CblasUpper && trans == CblasTrans) {
64 /* form B := alpha * (TriU(A))' *B */
66 for (i = n1; i > 0 && i--;) {
67 for (j = 0; j < n2; j++) {
70 for (k = 0; k < i; k++) {
71 temp += A[lda * k + i] * B[k * ldb + j];
75 temp += A[i * lda + i] * B[i * ldb + j];
77 temp += B[i * ldb + j];
80 B[ldb * i + j] = alpha * temp;
84 } else if (side == CblasLeft && uplo == CblasLower && trans == CblasNoTrans) {
86 /* form B := alpha * TriL(A)*B */
89 for (i = n1; i > 0 && i--;) {
90 for (j = 0; j < n2; j++) {
93 for (k = 0; k < i; k++) {
94 temp += A[lda * i + k] * B[k * ldb + j];
98 temp += A[i * lda + i] * B[i * ldb + j];
100 temp += B[i * ldb + j];
103 B[ldb * i + j] = alpha * temp;
109 } else if (side == CblasLeft && uplo == CblasLower && trans == CblasTrans) {
111 /* form B := alpha * TriL(A)' *B */
113 for (i = 0; i < n1; i++) {
114 for (j = 0; j < n2; j++) {
118 temp = A[i * lda + i] * B[i * ldb + j];
120 temp = B[i * ldb + j];
123 for (k = i + 1; k < n1; k++) {
124 temp += A[lda * k + i] * B[k * ldb + j];
127 B[ldb * i + j] = alpha * temp;
131 } else if (side == CblasRight && uplo == CblasUpper && trans == CblasNoTrans) {
133 /* form B := alpha * B * TriU(A) */
135 for (i = 0; i < n1; i++) {
136 for (j = n2; j > 0 && j--;) {
139 for (k = 0; k < j; k++) {
140 temp += A[lda * k + j] * B[i * ldb + k];
144 temp += A[j * lda + j] * B[i * ldb + j];
146 temp += B[i * ldb + j];
149 B[ldb * i + j] = alpha * temp;
153 } else if (side == CblasRight && uplo == CblasUpper && trans == CblasTrans) {
155 /* form B := alpha * B * (TriU(A))' */
157 for (i = 0; i < n1; i++) {
158 for (j = 0; j < n2; j++) {
162 temp = A[j * lda + j] * B[i * ldb + j];
164 temp = B[i * ldb + j];
167 for (k = j + 1; k < n2; k++) {
168 temp += A[lda * j + k] * B[i * ldb + k];
171 B[ldb * i + j] = alpha * temp;
175 } else if (side == CblasRight && uplo == CblasLower && trans == CblasNoTrans) {
177 /* form B := alpha *B * TriL(A) */
179 for (i = 0; i < n1; i++) {
180 for (j = 0; j < n2; j++) {
184 temp = A[j * lda + j] * B[i * ldb + j];
186 temp = B[i * ldb + j];
189 for (k = j + 1; k < n2; k++) {
190 temp += A[lda * k + j] * B[i * ldb + k];
194 B[ldb * i + j] = alpha * temp;
198 } else if (side == CblasRight && uplo == CblasLower && trans == CblasTrans) {
200 /* form B := alpha * B * TriL(A)' */
202 for (i = 0; i < n1; i++) {
203 for (j = n2; j > 0 && j--;) {
206 for (k = 0; k < j; k++) {
207 temp += A[lda * j + k] * B[i * ldb + k];
211 temp += A[j * lda + j] * B[i * ldb + j];
213 temp += B[i * ldb + j];
216 B[ldb * i + j] = alpha * temp;
221 BLAS_ERROR("unrecognized operation");