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