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

   1 /* blas/source_syrk_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   int uplo, trans;
  23
  24   const BASE alpha_real = CONST_REAL0(alpha);
  25   const BASE alpha_imag = CONST_IMAG0(alpha);
  26
  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     uplo = Uplo;
  36     /* FIXME: original blas does not make distinction between Trans and ConjTrans?? */
  37     trans = (Trans == CblasNoTrans) ? CblasNoTrans : CblasTrans;
  38   } else {
  39     uplo = (Uplo == CblasUpper) ? CblasLower : CblasUpper;
  40     trans = (Trans == CblasNoTrans) ? CblasTrans : CblasNoTrans;
  41   }
  42
  43   /* form  y := beta*y */
  44   if (beta_real == 0.0 && beta_imag == 0.0) {
  45     if (uplo == CblasUpper) {
  46       for (i = 0; i < N; i++) {
  47         for (j = i; j < N; j++) {
  48           REAL(C, ldc * i + j) = 0.0;
  49           IMAG(C, ldc * i + j) = 0.0;
  50         }
  51       }
  52     } else {
  53       for (i = 0; i < N; i++) {
  54         for (j = 0; j <= i; j++) {
  55           REAL(C, ldc * i + j) = 0.0;
  56           IMAG(C, ldc * i + j) = 0.0;
  57         }
  58       }
  59     }
  60   } else if (!(beta_real == 1.0 && beta_imag == 0.0)) {
  61     if (uplo == CblasUpper) {
  62       for (i = 0; i < N; i++) {
  63         for (j = i; j < N; j++) {
  64           const BASE Cij_real = REAL(C, ldc * i + j);
  65           const BASE Cij_imag = IMAG(C, ldc * i + j);
  66           REAL(C, ldc * i + j) = beta_real * Cij_real - beta_imag * Cij_imag;
  67           IMAG(C, ldc * i + j) = beta_real * Cij_imag + beta_imag * Cij_real;
  68         }
  69       }
  70     } else {
  71       for (i = 0; i < N; i++) {
  72         for (j = 0; j <= i; j++) {
  73           const BASE Cij_real = REAL(C, ldc * i + j);
  74           const BASE Cij_imag = IMAG(C, ldc * i + j);
  75           REAL(C, ldc * i + j) = beta_real * Cij_real - beta_imag * Cij_imag;
  76           IMAG(C, ldc * i + j) = beta_real * Cij_imag + beta_imag * Cij_real;
  77         }
  78       }
  79     }
  80   }
  81
  82   if (alpha_real == 0.0 && alpha_imag == 0.0)
  83     return;
  84
  85   if (uplo == CblasUpper && trans == CblasNoTrans) {
  86
  87     for (i = 0; i < N; i++) {
  88       for (j = i; j < N; j++) {
  89         BASE temp_real = 0.0;
  90         BASE temp_imag = 0.0;
  91         for (k = 0; k < K; k++) {
  92           const BASE Aik_real = CONST_REAL(A, i * lda + k);
  93           const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
  94           const BASE Ajk_real = CONST_REAL(A, j * lda + k);
  95           const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
  96           temp_real += Aik_real * Ajk_real - Aik_imag * Ajk_imag;
  97           temp_imag += Aik_real * Ajk_imag + Aik_imag * Ajk_real;
  98         }
  99         REAL(C, i * ldc + j) += alpha_real * temp_real - alpha_imag * temp_imag;
 100         IMAG(C, i * ldc + j) += alpha_real * temp_imag + alpha_imag * temp_real;
 101       }
 102     }
 103
 104   } else if (uplo == CblasUpper && trans == CblasTrans) {
 105
 106     for (i = 0; i < N; i++) {
 107       for (j = i; j < N; j++) {
 108         BASE temp_real = 0.0;
 109         BASE temp_imag = 0.0;
 110         for (k = 0; k < K; k++) {
 111           const BASE Aki_real = CONST_REAL(A, k * lda + i);
 112           const BASE Aki_imag = CONST_IMAG(A, k * lda + i);
 113           const BASE Akj_real = CONST_REAL(A, k * lda + j);
 114           const BASE Akj_imag = CONST_IMAG(A, k * lda + j);
 115           temp_real += Aki_real * Akj_real - Aki_imag * Akj_imag;
 116           temp_imag += Aki_real * Akj_imag + Aki_imag * Akj_real;
 117         }
 118         REAL(C, i * ldc + j) += alpha_real * temp_real - alpha_imag * temp_imag;
 119         IMAG(C, i * ldc + j) += alpha_real * temp_imag + alpha_imag * temp_real;
 120       }
 121     }
 122
 123   } else if (uplo == CblasLower && trans == CblasNoTrans) {
 124
 125     for (i = 0; i < N; i++) {
 126       for (j = 0; j <= i; j++) {
 127         BASE temp_real = 0.0;
 128         BASE temp_imag = 0.0;
 129         for (k = 0; k < K; k++) {
 130           const BASE Aik_real = CONST_REAL(A, i * lda + k);
 131           const BASE Aik_imag = CONST_IMAG(A, i * lda + k);
 132           const BASE Ajk_real = CONST_REAL(A, j * lda + k);
 133           const BASE Ajk_imag = CONST_IMAG(A, j * lda + k);
 134           temp_real += Aik_real * Ajk_real - Aik_imag * Ajk_imag;
 135           temp_imag += Aik_real * Ajk_imag + Aik_imag * Ajk_real;
 136         }
 137         REAL(C, i * ldc + j) += alpha_real * temp_real - alpha_imag * temp_imag;
 138         IMAG(C, i * ldc + j) += alpha_real * temp_imag + alpha_imag * temp_real;
 139       }
 140     }
 141
 142   } else if (uplo == CblasLower && trans == CblasTrans) {
 143
 144     for (i = 0; i < N; i++) {
 145       for (j = 0; j <= i; j++) {
 146         BASE temp_real = 0.0;
 147         BASE temp_imag = 0.0;
 148         for (k = 0; k < K; k++) {
 149           const BASE Aki_real = CONST_REAL(A, k * lda + i);
 150           const BASE Aki_imag = CONST_IMAG(A, k * lda + i);
 151           const BASE Akj_real = CONST_REAL(A, k * lda + j);
 152           const BASE Akj_imag = CONST_IMAG(A, k * lda + j);
 153           temp_real += Aki_real * Akj_real - Aki_imag * Akj_imag;
 154           temp_imag += Aki_real * Akj_imag + Aki_imag * Akj_real;
 155         }
 156         REAL(C, i * ldc + j) += alpha_real * temp_real - alpha_imag * temp_imag;
 157         IMAG(C, i * ldc + j) += alpha_real * temp_imag + alpha_imag * temp_real;
 158       }
 159     }
 160
 161   } else {
 162     BLAS_ERROR("unrecognized operation");
 163   }
 164 }