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

   1 /* blas/source_hemm.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   int uplo, side;
  24
  25   const BASE alpha_real = CONST_REAL0(alpha);
  26   const BASE alpha_imag = CONST_IMAG0(alpha);
  27
  28   const BASE beta_real = CONST_REAL0(beta);
  29   const BASE beta_imag = CONST_IMAG0(beta);
  30
  31   if ((alpha_real == 0.0 && alpha_imag == 0.0)
  32       && (beta_real == 1.0 && beta_imag == 0.0))
  33     return;
  34
  35   if (Order == CblasRowMajor) {
  36     n1 = M;
  37     n2 = N;
  38     uplo = Uplo;
  39     side = Side;
  40   } else {
  41     n1 = N;
  42     n2 = M;
  43     uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;
  44     side = (Side == CblasLeft) ? CblasRight : CblasLeft;
  45   }
  46
  47   /* form  y := beta*y */
  48   if (beta_real == 0.0 && beta_imag == 0.0) {
  49     for (i = 0; i < n1; i++) {
  50       for (j = 0; j < n2; j++) {
  51         REAL(C, ldc * i + j) = 0.0;
  52         IMAG(C, ldc * i + j) = 0.0;
  53       }
  54     }
  55   } else if (!(beta_real == 1.0 && beta_imag == 0.0)) {
  56     for (i = 0; i < n1; i++) {
  57       for (j = 0; j < n2; j++) {
  58         const BASE Cij_real = REAL(C, ldc * i + j);
  59         const BASE Cij_imag = IMAG(C, ldc * i + j);
  60         REAL(C, ldc * i + j) = beta_real * Cij_real - beta_imag * Cij_imag;
  61         IMAG(C, ldc * i + j) = beta_real * Cij_imag + beta_imag * Cij_real;
  62       }
  63     }
  64   }
  65
  66   if (alpha_real == 0.0 && alpha_imag == 0.0)
  67     return;
  68
  69   if (side == CblasLeft && uplo == CblasUpper) {
  70
  71     /* form  C := alpha*A*B + C */
  72
  73     for (i = 0; i < n1; i++) {
  74       for (j = 0; j < n2; j++) {
  75         const BASE Bij_real = CONST_REAL(B, ldb * i + j);
  76         const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
  77         const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
  78         const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
  79         BASE temp2_real = 0.0;
  80         BASE temp2_imag = 0.0;
  81         {
  82           const BASE Aii_real = CONST_REAL(A, i * lda + i);
  83           /* const BASE Aii_imag = 0.0; */
  84           REAL(C, i * ldc + j) += temp1_real * Aii_real;
  85           IMAG(C, i * ldc + j) += temp1_imag * Aii_real;
  86         }
  87         for (k = i + 1; k < n1; k++) {
  88           const BASE Aik_real = CONST_REAL(A, i * lda + k);
  89           const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
  90           const BASE Bkj_real = CONST_REAL(B, ldb * k + j);
  91           const BASE Bkj_imag = CONST_IMAG(B, ldb * k + j);
  92           REAL(C, k * ldc + j) += Aik_real * temp1_real - (-Aik_imag) * temp1_imag;
  93           IMAG(C, k * ldc + j) += Aik_real * temp1_imag + (-Aik_imag) * temp1_real;
  94           temp2_real += Aik_real * Bkj_real - Aik_imag * Bkj_imag;
  95           temp2_imag += Aik_real * Bkj_imag + Aik_imag * Bkj_real;
  96         }
  97         REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
  98         IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
  99       }
 100     }
 101
 102   } else if (side == CblasLeft && uplo == CblasLower) {
 103
 104     /* form  C := alpha*A*B + C */
 105
 106     for (i = 0; i < n1; i++) {
 107       for (j = 0; j < n2; j++) {
 108         const BASE Bij_real = CONST_REAL(B, ldb * i + j);
 109         const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
 110         const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
 111         const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
 112         BASE temp2_real = 0.0;
 113         BASE temp2_imag = 0.0;
 114         for (k = 0; k < i; k++) {
 115           const BASE Aik_real = CONST_REAL(A, i * lda + k);
 116           const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
 117           const BASE Bkj_real = CONST_REAL(B, ldb * k + j);
 118           const BASE Bkj_imag = CONST_IMAG(B, ldb * k + j);
 119           REAL(C, k * ldc + j) += Aik_real * temp1_real - (-Aik_imag) * temp1_imag;
 120           IMAG(C, k * ldc + j) += Aik_real * temp1_imag + (-Aik_imag) * temp1_real;
 121           temp2_real += Aik_real * Bkj_real - Aik_imag * Bkj_imag;
 122           temp2_imag += Aik_real * Bkj_imag + Aik_imag * Bkj_real;
 123         }
 124         {
 125           const BASE Aii_real = CONST_REAL(A, i * lda + i);
 126           /* const BASE Aii_imag = 0.0; */
 127           REAL(C, i * ldc + j) += temp1_real * Aii_real;
 128           IMAG(C, i * ldc + j) += temp1_imag * Aii_real;
 129         }
 130         REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
 131         IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
 132       }
 133     }
 134
 135   } else if (side == CblasRight && uplo == CblasUpper) {
 136
 137     /* form  C := alpha*B*A + C */
 138
 139     for (i = 0; i < n1; i++) {
 140       for (j = 0; j < n2; j++) {
 141         const BASE Bij_real = CONST_REAL(B, ldb * i + j);
 142         const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
 143         const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
 144         const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
 145         BASE temp2_real = 0.0;
 146         BASE temp2_imag = 0.0;
 147         {
 148           const BASE Ajj_real = CONST_REAL(A, j * lda + j);
 149           /* const BASE Ajj_imag = 0.0; */
 150           REAL(C, i * ldc + j) += temp1_real * Ajj_real;
 151           IMAG(C, i * ldc + j) += temp1_imag * Ajj_real;
 152         }
 153         for (k = j + 1; k < n2; k++) {
 154           const BASE Ajk_real = CONST_REAL(A, j * lda + k);
 155           const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
 156           const BASE Bik_real = CONST_REAL(B, ldb * i + k);
 157           const BASE Bik_imag = CONST_IMAG(B, ldb * i + k);
 158           REAL(C, i * ldc + k) += temp1_real * Ajk_real - temp1_imag * Ajk_imag;
 159           IMAG(C, i * ldc + k) += temp1_real * Ajk_imag + temp1_imag * Ajk_real;
 160           temp2_real += Bik_real * Ajk_real - Bik_imag * (-Ajk_imag);
 161           temp2_imag += Bik_real * (-Ajk_imag) + Bik_imag * Ajk_real;
 162         }
 163         REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
 164         IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
 165       }
 166     }
 167
 168   } else if (side == CblasRight && uplo == CblasLower) {
 169
 170     /* form  C := alpha*B*A + C */
 171
 172     for (i = 0; i < n1; i++) {
 173       for (j = 0; j < n2; j++) {
 174         const BASE Bij_real = CONST_REAL(B, ldb * i + j);
 175         const BASE Bij_imag = CONST_IMAG(B, ldb * i + j);
 176         const BASE temp1_real = alpha_real * Bij_real - alpha_imag * Bij_imag;
 177         const BASE temp1_imag = alpha_real * Bij_imag + alpha_imag * Bij_real;
 178         BASE temp2_real = 0.0;
 179         BASE temp2_imag = 0.0;
 180         for (k = 0; k < j; k++) {
 181           const BASE Ajk_real = CONST_REAL(A, j * lda + k);
 182           const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
 183           const BASE Bik_real = CONST_REAL(B, ldb * i + k);
 184           const BASE Bik_imag = CONST_IMAG(B, ldb * i + k);
 185           REAL(C, i * ldc + k) += temp1_real * Ajk_real - temp1_imag * Ajk_imag;
 186           IMAG(C, i * ldc + k) += temp1_real * Ajk_imag + temp1_imag * Ajk_real;
 187           temp2_real += Bik_real * Ajk_real - Bik_imag * (-Ajk_imag);
 188           temp2_imag += Bik_real * (-Ajk_imag) + Bik_imag * Ajk_real;
 189         }
 190         {
 191           const BASE Ajj_real = CONST_REAL(A, j * lda + j);
 192           /* const BASE Ajj_imag = 0.0; */
 193           REAL(C, i * ldc + j) += temp1_real * Ajj_real;
 194           IMAG(C, i * ldc + j) += temp1_imag * Ajj_real;
 195         }
 196         REAL(C, i * ldc + j) += alpha_real * temp2_real - alpha_imag * temp2_imag;
 197         IMAG(C, i * ldc + j) += alpha_real * temp2_imag + alpha_imag * temp2_real;
 198       }
 199     }
 200
 201   } else {
 202     BLAS_ERROR("unrecognized operation");
 203   }
 204 }