1 /* blas/source_trsm_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 * inv(TriU(A)) *B */ 43 44 if (alpha != 1.0) { 45 for (i = 0; i < n1; i++) { 46 for (j = 0; j < n2; j++) { 47 B[ldb * i + j] *= alpha; 48 } 49 } 50 } 51 52 for (i = n1; i > 0 && i--;) { 53 if (nonunit) { 54 BASE Aii = A[lda * i + i]; 55 for (j = 0; j < n2; j++) { 56 B[ldb * i + j] /= Aii; 57 } 58 } 59 60 for (k = 0; k < i; k++) { 61 const BASE Aki = A[k * lda + i]; 62 for (j = 0; j < n2; j++) { 63 B[ldb * k + j] -= Aki * B[ldb * i + j]; 64 } 65 } 66 } 67 68 } else if (side == CblasLeft && uplo == CblasUpper && trans == CblasTrans) { 69 70 /* form B := alpha * inv(TriU(A))' *B */ 71 72 if (alpha != 1.0) { 73 for (i = 0; i < n1; i++) { 74 for (j = 0; j < n2; j++) { 75 B[ldb * i + j] *= alpha; 76 } 77 } 78 } 79 80 for (i = 0; i < n1; i++) { 81 if (nonunit) { 82 BASE Aii = A[lda * i + i]; 83 for (j = 0; j < n2; j++) { 84 B[ldb * i + j] /= Aii; 85 } 86 } 87 88 for (k = i + 1; k < n1; k++) { 89 const BASE Aik = A[i * lda + k]; 90 for (j = 0; j < n2; j++) { 91 B[ldb * k + j] -= Aik * B[ldb * i + j]; 92 } 93 } 94 } 95 96 } else if (side == CblasLeft && uplo == CblasLower && trans == CblasNoTrans) { 97 98 /* form B := alpha * inv(TriL(A))*B */ 99 100 101 if (alpha != 1.0) { 102 for (i = 0; i < n1; i++) { 103 for (j = 0; j < n2; j++) { 104 B[ldb * i + j] *= alpha; 105 } 106 } 107 } 108 109 for (i = 0; i < n1; i++) { 110 if (nonunit) { 111 BASE Aii = A[lda * i + i]; 112 for (j = 0; j < n2; j++) { 113 B[ldb * i + j] /= Aii; 114 } 115 } 116 117 for (k = i + 1; k < n1; k++) { 118 const BASE Aki = A[k * lda + i]; 119 for (j = 0; j < n2; j++) { 120 B[ldb * k + j] -= Aki * B[ldb * i + j]; 121 } 122 } 123 } 124 125 126 } else if (side == CblasLeft && uplo == CblasLower && trans == CblasTrans) { 127 128 /* form B := alpha * TriL(A)' *B */ 129 130 if (alpha != 1.0) { 131 for (i = 0; i < n1; i++) { 132 for (j = 0; j < n2; j++) { 133 B[ldb * i + j] *= alpha; 134 } 135 } 136 } 137 138 for (i = n1; i > 0 && i--;) { 139 if (nonunit) { 140 BASE Aii = A[lda * i + i]; 141 for (j = 0; j < n2; j++) { 142 B[ldb * i + j] /= Aii; 143 } 144 } 145 146 for (k = 0; k < i; k++) { 147 const BASE Aik = A[i * lda + k]; 148 for (j = 0; j < n2; j++) { 149 B[ldb * k + j] -= Aik * B[ldb * i + j]; 150 } 151 } 152 } 153 154 } else if (side == CblasRight && uplo == CblasUpper && trans == CblasNoTrans) { 155 156 /* form B := alpha * B * inv(TriU(A)) */ 157 158 if (alpha != 1.0) { 159 for (i = 0; i < n1; i++) { 160 for (j = 0; j < n2; j++) { 161 B[ldb * i + j] *= alpha; 162 } 163 } 164 } 165 166 for (i = 0; i < n1; i++) { 167 for (j = 0; j < n2; j++) { 168 if (nonunit) { 169 BASE Ajj = A[lda * j + j]; 170 B[ldb * i + j] /= Ajj; 171 } 172 173 { 174 BASE Bij = B[ldb * i + j]; 175 for (k = j + 1; k < n2; k++) { 176 B[ldb * i + k] -= A[j * lda + k] * Bij; 177 } 178 } 179 } 180 } 181 182 } else if (side == CblasRight && uplo == CblasUpper && trans == CblasTrans) { 183 184 /* form B := alpha * B * inv(TriU(A))' */ 185 186 if (alpha != 1.0) { 187 for (i = 0; i < n1; i++) { 188 for (j = 0; j < n2; j++) { 189 B[ldb * i + j] *= alpha; 190 } 191 } 192 } 193 194 for (i = 0; i < n1; i++) { 195 for (j = n2; j > 0 && j--;) { 196 197 if (nonunit) { 198 BASE Ajj = A[lda * j + j]; 199 B[ldb * i + j] /= Ajj; 200 } 201 202 { 203 BASE Bij = B[ldb * i + j]; 204 for (k = 0; k < j; k++) { 205 B[ldb * i + k] -= A[k * lda + j] * Bij; 206 } 207 } 208 } 209 } 210 211 212 } else if (side == CblasRight && uplo == CblasLower && trans == CblasNoTrans) { 213 214 /* form B := alpha * B * inv(TriL(A)) */ 215 216 if (alpha != 1.0) { 217 for (i = 0; i < n1; i++) { 218 for (j = 0; j < n2; j++) { 219 B[ldb * i + j] *= alpha; 220 } 221 } 222 } 223 224 for (i = 0; i < n1; i++) { 225 for (j = n2; j > 0 && j--;) { 226 227 if (nonunit) { 228 BASE Ajj = A[lda * j + j]; 229 B[ldb * i + j] /= Ajj; 230 } 231 232 { 233 BASE Bij = B[ldb * i + j]; 234 for (k = 0; k < j; k++) { 235 B[ldb * i + k] -= A[j * lda + k] * Bij; 236 } 237 } 238 } 239 } 240 241 } else if (side == CblasRight && uplo == CblasLower && trans == CblasTrans) { 242 243 /* form B := alpha * B * inv(TriL(A))' */ 244 245 246 if (alpha != 1.0) { 247 for (i = 0; i < n1; i++) { 248 for (j = 0; j < n2; j++) { 249 B[ldb * i + j] *= alpha; 250 } 251 } 252 } 253 254 for (i = 0; i < n1; i++) { 255 for (j = 0; j < n2; j++) { 256 if (nonunit) { 257 BASE Ajj = A[lda * j + j]; 258 B[ldb * i + j] /= Ajj; 259 } 260 261 { 262 BASE Bij = B[ldb * i + j]; 263 for (k = j + 1; k < n2; k++) { 264 B[ldb * i + k] -= A[k * lda + j] * Bij; 265 } 266 } 267 } 268 } 269 270 271 272 } else { 273 BLAS_ERROR("unrecognized operation"); 274 } 275 } 276