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

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