1 /* blas/source_trmm_r.h
2  *
3  * Copyright (C) 2001, 2007 Brian Gough
4  *
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.
9  *
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.
14  *
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.
18  */
19 
20 {
21   INDEX i, j, k;
22   INDEX n1, n2;
23   const int nonunit = (Diag == CblasNonUnit);
24   int side, uplo, trans;
25 
26   if (Order == CblasRowMajor) {
27     n1 = M;
28     n2 = N;
29     side = Side;
30     uplo = Uplo;
31     trans = (TransA == CblasConjTrans) ? CblasTrans : TransA;
32   } else {
33     n1 = N;
34     n2 = M;
35     side = (Side == CblasLeft) ? CblasRight : CblasLeft;
36     uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;
37     trans = (TransA == CblasConjTrans) ? CblasTrans : TransA;
38   }
39 
40   if (side == CblasLeft && uplo == CblasUpper && trans == CblasNoTrans) {
41 
42     /* form  B := alpha * TriU(A)*B */
43 
44     for (i = 0; i < n1; i++) {
45       for (j = 0; j < n2; j++) {
46         BASE temp = 0.0;
47 
48         if (nonunit) {
49           temp = A[i * lda + i] * B[i * ldb + j];
50         } else {
51           temp = B[i * ldb + j];
52         }
53 
54         for (k = i + 1; k < n1; k++) {
55           temp += A[lda * i + k] * B[k * ldb + j];
56         }
57 
58         B[ldb * i + j] = alpha * temp;
59       }
60     }
61 
62   } else if (side == CblasLeft && uplo == CblasUpper && trans == CblasTrans) {
63 
64     /* form  B := alpha * (TriU(A))' *B */
65 
66     for (i = n1; i > 0 && i--;) {
67       for (j = 0; j < n2; j++) {
68         BASE temp = 0.0;
69 
70         for (k = 0; k < i; k++) {
71           temp += A[lda * k + i] * B[k * ldb + j];
72         }
73 
74         if (nonunit) {
75           temp += A[i * lda + i] * B[i * ldb + j];
76         } else {
77           temp += B[i * ldb + j];
78         }
79 
80         B[ldb * i + j] = alpha * temp;
81       }
82     }
83 
84   } else if (side == CblasLeft && uplo == CblasLower && trans == CblasNoTrans) {
85 
86     /* form  B := alpha * TriL(A)*B */
87 
88 
89     for (i = n1; i > 0 && i--;) {
90       for (j = 0; j < n2; j++) {
91         BASE temp = 0.0;
92 
93         for (k = 0; k < i; k++) {
94           temp += A[lda * i + k] * B[k * ldb + j];
95         }
96 
97         if (nonunit) {
98           temp += A[i * lda + i] * B[i * ldb + j];
99         } else {
100           temp += B[i * ldb + j];
101         }
102 
103         B[ldb * i + j] = alpha * temp;
104       }
105     }
106 
107 
108 
109   } else if (side == CblasLeft && uplo == CblasLower && trans == CblasTrans) {
110 
111     /* form  B := alpha * TriL(A)' *B */
112 
113     for (i = 0; i < n1; i++) {
114       for (j = 0; j < n2; j++) {
115         BASE temp = 0.0;
116 
117         if (nonunit) {
118           temp = A[i * lda + i] * B[i * ldb + j];
119         } else {
120           temp = B[i * ldb + j];
121         }
122 
123         for (k = i + 1; k < n1; k++) {
124           temp += A[lda * k + i] * B[k * ldb + j];
125         }
126 
127         B[ldb * i + j] = alpha * temp;
128       }
129     }
130 
131   } else if (side == CblasRight && uplo == CblasUpper && trans == CblasNoTrans) {
132 
133     /* form  B := alpha * B * TriU(A) */
134 
135     for (i = 0; i < n1; i++) {
136       for (j = n2; j > 0 && j--;) {
137         BASE temp = 0.0;
138 
139         for (k = 0; k < j; k++) {
140           temp += A[lda * k + j] * B[i * ldb + k];
141         }
142 
143         if (nonunit) {
144           temp += A[j * lda + j] * B[i * ldb + j];
145         } else {
146           temp += B[i * ldb + j];
147         }
148 
149         B[ldb * i + j] = alpha * temp;
150       }
151     }
152 
153   } else if (side == CblasRight && uplo == CblasUpper && trans == CblasTrans) {
154 
155     /* form  B := alpha * B * (TriU(A))' */
156 
157     for (i = 0; i < n1; i++) {
158       for (j = 0; j < n2; j++) {
159         BASE temp = 0.0;
160 
161         if (nonunit) {
162           temp = A[j * lda + j] * B[i * ldb + j];
163         } else {
164           temp = B[i * ldb + j];
165         }
166 
167         for (k = j + 1; k < n2; k++) {
168           temp += A[lda * j + k] * B[i * ldb + k];
169         }
170 
171         B[ldb * i + j] = alpha * temp;
172       }
173     }
174 
175   } else if (side == CblasRight && uplo == CblasLower && trans == CblasNoTrans) {
176 
177     /* form  B := alpha *B * TriL(A) */
178 
179     for (i = 0; i < n1; i++) {
180       for (j = 0; j < n2; j++) {
181         BASE temp = 0.0;
182 
183         if (nonunit) {
184           temp = A[j * lda + j] * B[i * ldb + j];
185         } else {
186           temp = B[i * ldb + j];
187         }
188 
189         for (k = j + 1; k < n2; k++) {
190           temp += A[lda * k + j] * B[i * ldb + k];
191         }
192 
193 
194         B[ldb * i + j] = alpha * temp;
195       }
196     }
197 
198   } else if (side == CblasRight && uplo == CblasLower && trans == CblasTrans) {
199 
200     /* form  B := alpha * B * TriL(A)' */
201 
202     for (i = 0; i < n1; i++) {
203       for (j = n2; j > 0 && j--;) {
204         BASE temp = 0.0;
205 
206         for (k = 0; k < j; k++) {
207           temp += A[lda * j + k] * B[i * ldb + k];
208         }
209 
210         if (nonunit) {
211           temp += A[j * lda + j] * B[i * ldb + j];
212         } else {
213           temp += B[i * ldb + j];
214         }
215 
216         B[ldb * i + j] = alpha * temp;
217       }
218     }
219 
220   } else {
221     BLAS_ERROR("unrecognized operation");
222   }
223 }
224