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

   1 /* blas/source_tbmv_c.h
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   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   const int conj = (TransA == CblasConjTrans) ? -1 : 1;
  22   const int Trans = (TransA != CblasConjTrans) ? TransA : CblasTrans;
  23   const int nonunit = (Diag == CblasNonUnit);
  24   INDEX i, j;
  25
  26
  27   if ((order == CblasRowMajor && Trans == CblasNoTrans && Uplo == CblasUpper)
  28       || (order == CblasColMajor && Trans == CblasTrans && Uplo == CblasLower)) {
  29     /* form  x := A*x */
  30
  31     INDEX ix = OFFSET(N, incX);
  32     for (i = 0; i < N; i++) {
  33       BASE temp_r = 0.0;
  34       BASE temp_i = 0.0;
  35       const INDEX j_min = i + 1;
  36       const INDEX j_max = GSL_MIN(N, i + K + 1);
  37       INDEX jx = OFFSET(N, incX) + incX * j_min;
  38       for (j = j_min; j < j_max; j++) {
  39         const BASE x_real = REAL(X, jx);
  40         const BASE x_imag = IMAG(X, jx);
  41         const BASE A_real = CONST_REAL(A, lda * i + (j - i));
  42         const BASE A_imag = conj * CONST_IMAG(A, lda * i + (j - i));
  43
  44         temp_r += A_real * x_real - A_imag * x_imag;
  45         temp_i += A_real * x_imag + A_imag * x_real;
  46
  47         jx += incX;
  48       }
  49       if (nonunit) {
  50         const BASE x_real = REAL(X, ix);
  51         const BASE x_imag = IMAG(X, ix);
  52         const BASE A_real = CONST_REAL(A, lda * i + 0);
  53         const BASE A_imag = conj * CONST_IMAG(A, lda * i + 0);
  54
  55         REAL(X, ix) = temp_r + (A_real * x_real - A_imag * x_imag);
  56         IMAG(X, ix) = temp_i + (A_real * x_imag + A_imag * x_real);
  57       } else {
  58         REAL(X, ix) += temp_r;
  59         IMAG(X, ix) += temp_i;
  60       }
  61       ix += incX;
  62     }
  63   } else if ((order == CblasRowMajor && Trans == CblasNoTrans && Uplo == CblasLower)
  64              || (order == CblasColMajor && Trans == CblasTrans && Uplo == CblasUpper)) {
  65     INDEX ix = OFFSET(N, incX) + (N - 1) * incX;
  66
  67     for (i = N; i > 0 && i--;) {        /*  N-1 ... 0 */
  68       BASE temp_r = 0.0;
  69       BASE temp_i = 0.0;
  70       const INDEX j_min = (K > i ? 0 : i - K);
  71       const INDEX j_max = i;
  72       INDEX jx = OFFSET(N, incX) + j_min * incX;
  73       for (j = j_min; j < j_max; j++) {
  74         const BASE x_real = REAL(X, jx);
  75         const BASE x_imag = IMAG(X, jx);
  76         const BASE A_real = CONST_REAL(A, lda * i + (K - i + j));
  77         const BASE A_imag = conj * CONST_IMAG(A, lda * i + (K - i + j));
  78
  79         temp_r += A_real * x_real - A_imag * x_imag;
  80         temp_i += A_real * x_imag + A_imag * x_real;
  81
  82         jx += incX;
  83       }
  84       if (nonunit) {
  85         const BASE x_real = REAL(X, ix);
  86         const BASE x_imag = IMAG(X, ix);
  87         const BASE A_real = CONST_REAL(A, lda * i + K);
  88         const BASE A_imag = conj * CONST_IMAG(A, lda * i + K);
  89
  90         REAL(X, ix) = temp_r + (A_real * x_real - A_imag * x_imag);
  91         IMAG(X, ix) = temp_i + (A_real * x_imag + A_imag * x_real);
  92       } else {
  93         REAL(X, ix) += temp_r;
  94         IMAG(X, ix) += temp_i;
  95       }
  96       ix -= incX;
  97     }
  98   } else if ((order == CblasRowMajor && Trans == CblasTrans && Uplo == CblasUpper)
  99              || (order == CblasColMajor && Trans == CblasNoTrans && Uplo == CblasLower)) {
 100     /* form  x := A'*x */
 101
 102     INDEX ix = OFFSET(N, incX) + (N - 1) * incX;
 103     for (i = N; i > 0 && i--;) {        /*  N-1 ... 0 */
 104       BASE temp_r = 0.0;
 105       BASE temp_i = 0.0;
 106       const INDEX j_min = (K > i ? 0 : i - K);
 107       const INDEX j_max = i;
 108       INDEX jx = OFFSET(N, incX) + j_min * incX;
 109       for (j = j_min; j < j_max; j++) {
 110         const BASE x_real = REAL(X, jx);
 111         const BASE x_imag = IMAG(X, jx);
 112         const BASE A_real = CONST_REAL(A, lda * j + (i - j));
 113         const BASE A_imag = conj * CONST_IMAG(A, lda * j + (i - j));
 114
 115         temp_r += A_real * x_real - A_imag * x_imag;
 116         temp_i += A_real * x_imag + A_imag * x_real;
 117
 118         jx += incX;
 119       }
 120       if (nonunit) {
 121         const BASE x_real = REAL(X, ix);
 122         const BASE x_imag = IMAG(X, ix);
 123         const BASE A_real = CONST_REAL(A, lda * i + 0);
 124         const BASE A_imag = conj * CONST_IMAG(A, lda * i + 0);
 125
 126         REAL(X, ix) = temp_r + (A_real * x_real - A_imag * x_imag);
 127         IMAG(X, ix) = temp_i + (A_real * x_imag + A_imag * x_real);
 128       } else {
 129         REAL(X, ix) += temp_r;
 130         IMAG(X, ix) += temp_i;
 131       }
 132       ix -= incX;
 133     }
 134   } else if ((order == CblasRowMajor && Trans == CblasTrans && Uplo == CblasLower)
 135              || (order == CblasColMajor && Trans == CblasNoTrans && Uplo == CblasUpper)) {
 136     INDEX ix = OFFSET(N, incX);
 137     for (i = 0; i < N; i++) {
 138       BASE temp_r = 0.0;
 139       BASE temp_i = 0.0;
 140       const INDEX j_min = i + 1;
 141       const INDEX j_max = GSL_MIN(N, i + K + 1);
 142       INDEX jx = OFFSET(N, incX) + j_min * incX;
 143       for (j = j_min; j < j_max; j++) {
 144         const BASE x_real = REAL(X, jx);
 145         const BASE x_imag = IMAG(X, jx);
 146         const BASE A_real = CONST_REAL(A, lda * j + (K - j + i));
 147         const BASE A_imag = conj * CONST_IMAG(A, lda * j + (K - j + i));
 148
 149         temp_r += A_real * x_real - A_imag * x_imag;
 150         temp_i += A_real * x_imag + A_imag * x_real;
 151
 152         jx += incX;
 153       }
 154       if (nonunit) {
 155         const BASE x_real = REAL(X, ix);
 156         const BASE x_imag = IMAG(X, ix);
 157         const BASE A_real = CONST_REAL(A, lda * i + K);
 158         const BASE A_imag = conj * CONST_IMAG(A, lda * i + K);
 159
 160         REAL(X, ix) = temp_r + (A_real * x_real - A_imag * x_imag);
 161         IMAG(X, ix) = temp_i + (A_real * x_imag + A_imag * x_real);
 162       } else {
 163         REAL(X, ix) += temp_r;
 164         IMAG(X, ix) += temp_i;
 165       }
 166       ix += incX;
 167     }
 168   } else {
 169     BLAS_ERROR("unrecognized operation");
 170   }
 171 }