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