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

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