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

   1 /* blas/source_trmm_c.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   const int conj = (TransA == CblasConjTrans) ? -1 : 1;
  25   int side, uplo, trans;
  26
  27   const BASE alpha_real = CONST_REAL0(alpha);
  28   const BASE alpha_imag = CONST_IMAG0(alpha);
  29
  30   if (Order == CblasRowMajor) {
  31     n1 = M;
  32     n2 = N;
  33     side = Side;
  34     uplo = Uplo;
  35     trans = (TransA == CblasNoTrans) ? CblasNoTrans : CblasTrans;
  36   } else {
  37     n1 = N;
  38     n2 = M;
  39     side = (Side == CblasLeft) ? CblasRight : CblasLeft;        /* exchanged */
  40     uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;      /* exchanged */
  41     trans = (TransA == CblasNoTrans) ? CblasNoTrans : CblasTrans;       /* same */
  42   }
  43
  44   if (side == CblasLeft && uplo == CblasUpper && trans == CblasNoTrans) {
  45
  46     /* form  B := alpha * TriU(A)*B */
  47
  48     for (i = 0; i < n1; i++) {
  49       for (j = 0; j < n2; j++) {
  50         BASE temp_real = 0.0;
  51         BASE temp_imag = 0.0;
  52
  53         if (nonunit) {
  54           const BASE Aii_real = CONST_REAL(A, i * lda + i);
  55           const BASE Aii_imag = conj * CONST_IMAG(A, i * lda + i);
  56           const BASE Bij_real = REAL(B, i * ldb + j);
  57           const BASE Bij_imag = IMAG(B, i * ldb + j);
  58           temp_real = Aii_real * Bij_real - Aii_imag * Bij_imag;
  59           temp_imag = Aii_real * Bij_imag + Aii_imag * Bij_real;
  60         } else {
  61           temp_real = REAL(B, i * ldb + j);
  62           temp_imag = IMAG(B, i * ldb + j);
  63         }
  64
  65         for (k = i + 1; k < n1; k++) {
  66           const BASE Aik_real = CONST_REAL(A, i * lda + k);
  67           const BASE Aik_imag = conj * CONST_IMAG(A, i * lda + k);
  68           const BASE Bkj_real = REAL(B, k * ldb + j);
  69           const BASE Bkj_imag = IMAG(B, k * ldb + j);
  70           temp_real += Aik_real * Bkj_real - Aik_imag * Bkj_imag;
  71           temp_imag += Aik_real * Bkj_imag + Aik_imag * Bkj_real;
  72         }
  73
  74         REAL(B, ldb * i + j) = alpha_real * temp_real - alpha_imag * temp_imag;
  75         IMAG(B, ldb * i + j) = alpha_real * temp_imag + alpha_imag * temp_real;
  76       }
  77     }
  78
  79   } else if (side == CblasLeft && uplo == CblasUpper && trans == CblasTrans) {
  80
  81     /* form  B := alpha * (TriU(A))' *B */
  82
  83     for (i = n1; i > 0 && i--;) {
  84       for (j = 0; j < n2; j++) {
  85         BASE temp_real = 0.0;
  86         BASE temp_imag = 0.0;
  87
  88         for (k = 0; k < i; k++) {
  89           const BASE Aki_real = CONST_REAL(A, k * lda + i);
  90           const BASE Aki_imag = conj * CONST_IMAG(A, k * lda + i);
  91           const BASE Bkj_real = REAL(B, k * ldb + j);
  92           const BASE Bkj_imag = IMAG(B, k * ldb + j);
  93           temp_real += Aki_real * Bkj_real - Aki_imag * Bkj_imag;
  94           temp_imag += Aki_real * Bkj_imag + Aki_imag * Bkj_real;
  95         }
  96
  97         if (nonunit) {
  98           const BASE Aii_real = CONST_REAL(A, i * lda + i);
  99           const BASE Aii_imag = conj * CONST_IMAG(A, i * lda + i);
 100           const BASE Bij_real = REAL(B, i * ldb + j);
 101           const BASE Bij_imag = IMAG(B, i * ldb + j);
 102           temp_real += Aii_real * Bij_real - Aii_imag * Bij_imag;
 103           temp_imag += Aii_real * Bij_imag + Aii_imag * Bij_real;
 104         } else {
 105           temp_real += REAL(B, i * ldb + j);
 106           temp_imag += IMAG(B, i * ldb + j);
 107         }
 108
 109         REAL(B, ldb * i + j) = alpha_real * temp_real - alpha_imag * temp_imag;
 110         IMAG(B, ldb * i + j) = alpha_real * temp_imag + alpha_imag * temp_real;
 111       }
 112     }
 113
 114   } else if (side == CblasLeft && uplo == CblasLower && trans == CblasNoTrans) {
 115
 116     /* form  B := alpha * TriL(A)*B */
 117
 118
 119     for (i = n1; i > 0 && i--;) {
 120       for (j = 0; j < n2; j++) {
 121         BASE temp_real = 0.0;
 122         BASE temp_imag = 0.0;
 123
 124         for (k = 0; k < i; k++) {
 125           const BASE Aik_real = CONST_REAL(A, i * lda + k);
 126           const BASE Aik_imag = conj * CONST_IMAG(A, i * lda + k);
 127           const BASE Bkj_real = REAL(B, k * ldb + j);
 128           const BASE Bkj_imag = IMAG(B, k * ldb + j);
 129           temp_real += Aik_real * Bkj_real - Aik_imag * Bkj_imag;
 130           temp_imag += Aik_real * Bkj_imag + Aik_imag * Bkj_real;
 131         }
 132
 133         if (nonunit) {
 134           const BASE Aii_real = CONST_REAL(A, i * lda + i);
 135           const BASE Aii_imag = conj * CONST_IMAG(A, i * lda + i);
 136           const BASE Bij_real = REAL(B, i * ldb + j);
 137           const BASE Bij_imag = IMAG(B, i * ldb + j);
 138           temp_real += Aii_real * Bij_real - Aii_imag * Bij_imag;
 139           temp_imag += Aii_real * Bij_imag + Aii_imag * Bij_real;
 140         } else {
 141           temp_real += REAL(B, i * ldb + j);
 142           temp_imag += IMAG(B, i * ldb + j);
 143         }
 144
 145         REAL(B, ldb * i + j) = alpha_real * temp_real - alpha_imag * temp_imag;
 146         IMAG(B, ldb * i + j) = alpha_real * temp_imag + alpha_imag * temp_real;
 147       }
 148     }
 149
 150
 151
 152   } else if (side == CblasLeft && uplo == CblasLower && trans == CblasTrans) {
 153
 154     /* form  B := alpha * TriL(A)' *B */
 155
 156     for (i = 0; i < n1; i++) {
 157       for (j = 0; j < n2; j++) {
 158         BASE temp_real = 0.0;
 159         BASE temp_imag = 0.0;
 160
 161         if (nonunit) {
 162           const BASE Aii_real = CONST_REAL(A, i * lda + i);
 163           const BASE Aii_imag = conj * CONST_IMAG(A, i * lda + i);
 164           const BASE Bij_real = REAL(B, i * ldb + j);
 165           const BASE Bij_imag = IMAG(B, i * ldb + j);
 166           temp_real = Aii_real * Bij_real - Aii_imag * Bij_imag;
 167           temp_imag = Aii_real * Bij_imag + Aii_imag * Bij_real;
 168         } else {
 169           temp_real = REAL(B, i * ldb + j);
 170           temp_imag = IMAG(B, i * ldb + j);
 171         }
 172
 173         for (k = i + 1; k < n1; k++) {
 174           const BASE Aki_real = CONST_REAL(A, k * lda + i);
 175           const BASE Aki_imag = conj * CONST_IMAG(A, k * lda + i);
 176           const BASE Bkj_real = REAL(B, k * ldb + j);
 177           const BASE Bkj_imag = IMAG(B, k * ldb + j);
 178           temp_real += Aki_real * Bkj_real - Aki_imag * Bkj_imag;
 179           temp_imag += Aki_real * Bkj_imag + Aki_imag * Bkj_real;
 180         }
 181
 182         REAL(B, ldb * i + j) = alpha_real * temp_real - alpha_imag * temp_imag;
 183         IMAG(B, ldb * i + j) = alpha_real * temp_imag + alpha_imag * temp_real;
 184       }
 185     }
 186
 187   } else if (side == CblasRight && uplo == CblasUpper && trans == CblasNoTrans) {
 188
 189     /* form  B := alpha * B * TriU(A) */
 190
 191     for (i = 0; i < n1; i++) {
 192       for (j = n2; j > 0 && j--;) {
 193         BASE temp_real = 0.0;
 194         BASE temp_imag = 0.0;
 195
 196         for (k = 0; k < j; k++) {
 197           const BASE Akj_real = CONST_REAL(A, k * lda + j);
 198           const BASE Akj_imag = conj * CONST_IMAG(A, k * lda + j);
 199           const BASE Bik_real = REAL(B, i * ldb + k);
 200           const BASE Bik_imag = IMAG(B, i * ldb + k);
 201           temp_real += Akj_real * Bik_real - Akj_imag * Bik_imag;
 202           temp_imag += Akj_real * Bik_imag + Akj_imag * Bik_real;
 203         }
 204
 205         if (nonunit) {
 206           const BASE Ajj_real = CONST_REAL(A, j * lda + j);
 207           const BASE Ajj_imag = conj * CONST_IMAG(A, j * lda + j);
 208           const BASE Bij_real = REAL(B, i * ldb + j);
 209           const BASE Bij_imag = IMAG(B, i * ldb + j);
 210           temp_real += Ajj_real * Bij_real - Ajj_imag * Bij_imag;
 211           temp_imag += Ajj_real * Bij_imag + Ajj_imag * Bij_real;
 212         } else {
 213           temp_real += REAL(B, i * ldb + j);
 214           temp_imag += IMAG(B, i * ldb + j);
 215         }
 216
 217         REAL(B, ldb * i + j) = alpha_real * temp_real - alpha_imag * temp_imag;
 218         IMAG(B, ldb * i + j) = alpha_real * temp_imag + alpha_imag * temp_real;
 219       }
 220     }
 221
 222   } else if (side == CblasRight && uplo == CblasUpper && trans == CblasTrans) {
 223
 224     /* form  B := alpha * B * (TriU(A))' */
 225
 226     for (i = 0; i < n1; i++) {
 227       for (j = 0; j < n2; j++) {
 228         BASE temp_real = 0.0;
 229         BASE temp_imag = 0.0;
 230
 231         if (nonunit) {
 232           const BASE Ajj_real = CONST_REAL(A, j * lda + j);
 233           const BASE Ajj_imag = conj * CONST_IMAG(A, j * lda + j);
 234           const BASE Bij_real = REAL(B, i * ldb + j);
 235           const BASE Bij_imag = IMAG(B, i * ldb + j);
 236           temp_real = Ajj_real * Bij_real - Ajj_imag * Bij_imag;
 237           temp_imag = Ajj_real * Bij_imag + Ajj_imag * Bij_real;
 238         } else {
 239           temp_real = REAL(B, i * ldb + j);
 240           temp_imag = IMAG(B, i * ldb + j);
 241         }
 242
 243         for (k = j + 1; k < n2; k++) {
 244           const BASE Ajk_real = CONST_REAL(A, j * lda + k);
 245           const BASE Ajk_imag = conj * CONST_IMAG(A, j * lda + k);
 246           const BASE Bik_real = REAL(B, i * ldb + k);
 247           const BASE Bik_imag = IMAG(B, i * ldb + k);
 248           temp_real += Ajk_real * Bik_real - Ajk_imag * Bik_imag;
 249           temp_imag += Ajk_real * Bik_imag + Ajk_imag * Bik_real;
 250         }
 251
 252         REAL(B, ldb * i + j) = alpha_real * temp_real - alpha_imag * temp_imag;
 253         IMAG(B, ldb * i + j) = alpha_real * temp_imag + alpha_imag * temp_real;
 254       }
 255     }
 256
 257   } else if (side == CblasRight && uplo == CblasLower && trans == CblasNoTrans) {
 258
 259     /* form  B := alpha *B * TriL(A) */
 260
 261     for (i = 0; i < n1; i++) {
 262       for (j = 0; j < n2; j++) {
 263         BASE temp_real = 0.0;
 264         BASE temp_imag = 0.0;
 265
 266         if (nonunit) {
 267           const BASE Ajj_real = CONST_REAL(A, j * lda + j);
 268           const BASE Ajj_imag = conj * CONST_IMAG(A, j * lda + j);
 269           const BASE Bij_real = REAL(B, i * ldb + j);
 270           const BASE Bij_imag = IMAG(B, i * ldb + j);
 271           temp_real = Ajj_real * Bij_real - Ajj_imag * Bij_imag;
 272           temp_imag = Ajj_real * Bij_imag + Ajj_imag * Bij_real;
 273         } else {
 274           temp_real = REAL(B, i * ldb + j);
 275           temp_imag = IMAG(B, i * ldb + j);
 276         }
 277
 278         for (k = j + 1; k < n2; k++) {
 279           const BASE Akj_real = CONST_REAL(A, k * lda + j);
 280           const BASE Akj_imag = conj * CONST_IMAG(A, k * lda + j);
 281           const BASE Bik_real = REAL(B, i * ldb + k);
 282           const BASE Bik_imag = IMAG(B, i * ldb + k);
 283           temp_real += Akj_real * Bik_real - Akj_imag * Bik_imag;
 284           temp_imag += Akj_real * Bik_imag + Akj_imag * Bik_real;
 285         }
 286
 287         REAL(B, ldb * i + j) = alpha_real * temp_real - alpha_imag * temp_imag;
 288         IMAG(B, ldb * i + j) = alpha_real * temp_imag + alpha_imag * temp_real;
 289       }
 290     }
 291
 292   } else if (side == CblasRight && uplo == CblasLower && trans == CblasTrans) {
 293
 294     /* form  B := alpha * B * TriL(A)' */
 295
 296     for (i = 0; i < n1; i++) {
 297       for (j = n2; j > 0 && j--;) {
 298         BASE temp_real = 0.0;
 299         BASE temp_imag = 0.0;
 300
 301         for (k = 0; k < j; k++) {
 302           const BASE Ajk_real = CONST_REAL(A, j * lda + k);
 303           const BASE Ajk_imag = conj * CONST_IMAG(A, j * lda + k);
 304           const BASE Bik_real = REAL(B, i * ldb + k);
 305           const BASE Bik_imag = IMAG(B, i * ldb + k);
 306           temp_real += Ajk_real * Bik_real - Ajk_imag * Bik_imag;
 307           temp_imag += Ajk_real * Bik_imag + Ajk_imag * Bik_real;
 308         }
 309
 310         if (nonunit) {
 311           const BASE Ajj_real = CONST_REAL(A, j * lda + j);
 312           const BASE Ajj_imag = conj * CONST_IMAG(A, j * lda + j);
 313           const BASE Bij_real = REAL(B, i * ldb + j);
 314           const BASE Bij_imag = IMAG(B, i * ldb + j);
 315           temp_real += Ajj_real * Bij_real - Ajj_imag * Bij_imag;
 316           temp_imag += Ajj_real * Bij_imag + Ajj_imag * Bij_real;
 317         } else {
 318           temp_real += REAL(B, i * ldb + j);
 319           temp_imag += IMAG(B, i * ldb + j);
 320         }
 321
 322         REAL(B, ldb * i + j) = alpha_real * temp_real - alpha_imag * temp_imag;
 323         IMAG(B, ldb * i + j) = alpha_real * temp_imag + alpha_imag * temp_real;
 324       }
 325     }
 326
 327   } else {
 328     BLAS_ERROR("unrecognized operation");
 329   }
 330 }