htswanalysis/MACS/lib/gsl/gsl-1.11/cblas/source_trsm_r.h

   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 }