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

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