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

   1 /* blas/source_trsv_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   INDEX ix, jx;
  26
  27   if (N == 0)
  28     return;
  29
  30   /* form  x := inv( A )*x */
  31
  32   if ((order == CblasRowMajor && Trans == CblasNoTrans && Uplo == CblasUpper)
  33       || (order == CblasColMajor && Trans == CblasTrans && Uplo == CblasLower)) {
  34
  35     ix = OFFSET(N, incX) + incX * (N - 1);
  36
  37     if (nonunit) {
  38       const BASE a_real = CONST_REAL(A, lda * (N - 1) + (N - 1));
  39       const BASE a_imag = conj * CONST_IMAG(A, lda * (N - 1) + (N - 1));
  40       const BASE x_real = REAL(X, ix);
  41       const BASE x_imag = IMAG(X, ix);
  42       const BASE s = xhypot(a_real, a_imag);
  43       const BASE b_real = a_real / s;
  44       const BASE b_imag = a_imag / s;
  45       REAL(X, ix) = (x_real * b_real + x_imag * b_imag) / s;
  46       IMAG(X, ix) = (x_imag * b_real - b_imag * x_real) / s;
  47     }
  48
  49     ix -= incX;
  50
  51     for (i = N - 1; i > 0 && i--;) {
  52       BASE tmp_real = REAL(X, ix);
  53       BASE tmp_imag = IMAG(X, ix);
  54       jx = ix + incX;
  55       for (j = i + 1; j < N; j++) {
  56         const BASE Aij_real = CONST_REAL(A, lda * i + j);
  57         const BASE Aij_imag = conj * CONST_IMAG(A, lda * i + j);
  58         const BASE x_real = REAL(X, jx);
  59         const BASE x_imag = IMAG(X, jx);
  60         tmp_real -= Aij_real * x_real - Aij_imag * x_imag;
  61         tmp_imag -= Aij_real * x_imag + Aij_imag * x_real;
  62         jx += incX;
  63       }
  64
  65       if (nonunit) {
  66         const BASE a_real = CONST_REAL(A, lda * i + i);
  67         const BASE a_imag = conj * CONST_IMAG(A, lda * i + i);
  68         const BASE s = xhypot(a_real, a_imag);
  69         const BASE b_real = a_real / s;
  70         const BASE b_imag = a_imag / s;
  71         REAL(X, ix) = (tmp_real * b_real + tmp_imag * b_imag) / s;
  72         IMAG(X, ix) = (tmp_imag * b_real - tmp_real * b_imag) / s;
  73       } else {
  74         REAL(X, ix) = tmp_real;
  75         IMAG(X, ix) = tmp_imag;
  76       }
  77       ix -= incX;
  78     }
  79
  80   } else if ((order == CblasRowMajor && Trans == CblasNoTrans && Uplo == CblasLower)
  81              || (order == CblasColMajor && Trans == CblasTrans && Uplo == CblasUpper)) {
  82     /* forward substitution */
  83
  84     ix = OFFSET(N, incX);
  85
  86     if (nonunit) {
  87       const BASE a_real = CONST_REAL(A, lda * 0 + 0);
  88       const BASE a_imag = conj * CONST_IMAG(A, lda * 0 + 0);
  89       const BASE x_real = REAL(X, ix);
  90       const BASE x_imag = IMAG(X, ix);
  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) = (x_real * b_real + x_imag * b_imag) / s;
  95       IMAG(X, ix) = (x_imag * b_real - b_imag * x_real) / s;
  96     }
  97
  98     ix += incX;
  99
 100     for (i = 1; i < N; i++) {
 101       BASE tmp_real = REAL(X, ix);
 102       BASE tmp_imag = IMAG(X, ix);
 103       jx = OFFSET(N, incX);
 104       for (j = 0; j < i; j++) {
 105         const BASE Aij_real = CONST_REAL(A, lda * i + j);
 106         const BASE Aij_imag = conj * CONST_IMAG(A, lda * i + j);
 107         const BASE x_real = REAL(X, jx);
 108         const BASE x_imag = IMAG(X, jx);
 109         tmp_real -= Aij_real * x_real - Aij_imag * x_imag;
 110         tmp_imag -= Aij_real * x_imag + Aij_imag * x_real;
 111         jx += incX;
 112       }
 113       if (nonunit) {
 114         const BASE a_real = CONST_REAL(A, lda * i + i);
 115         const BASE a_imag = conj * CONST_IMAG(A, lda * i + i);
 116         const BASE s = xhypot(a_real, a_imag);
 117         const BASE b_real = a_real / s;
 118         const BASE b_imag = a_imag / s;
 119         REAL(X, ix) = (tmp_real * b_real + tmp_imag * b_imag) / s;
 120         IMAG(X, ix) = (tmp_imag * b_real - tmp_real * b_imag) / s;
 121       } else {
 122         REAL(X, ix) = tmp_real;
 123         IMAG(X, ix) = tmp_imag;
 124       }
 125       ix += incX;
 126     }
 127   } else if ((order == CblasRowMajor && Trans == CblasTrans && Uplo == CblasUpper)
 128              || (order == CblasColMajor && Trans == CblasNoTrans && Uplo == CblasLower)) {
 129     /* form  x := inv( A' )*x */
 130
 131     /* forward substitution */
 132
 133     ix = OFFSET(N, incX);
 134
 135     if (nonunit) {
 136       const BASE a_real = CONST_REAL(A, lda * 0 + 0);
 137       const BASE a_imag = conj * CONST_IMAG(A, lda * 0 + 0);
 138       const BASE x_real = REAL(X, ix);
 139       const BASE x_imag = IMAG(X, ix);
 140       const BASE s = xhypot(a_real, a_imag);
 141       const BASE b_real = a_real / s;
 142       const BASE b_imag = a_imag / s;
 143       REAL(X, ix) = (x_real * b_real + x_imag * b_imag) / s;
 144       IMAG(X, ix) = (x_imag * b_real - b_imag * x_real) / s;
 145     }
 146
 147     ix += incX;
 148
 149     for (i = 1; i < N; i++) {
 150       BASE tmp_real = REAL(X, ix);
 151       BASE tmp_imag = IMAG(X, ix);
 152       jx = OFFSET(N, incX);
 153       for (j = 0; j < i; j++) {
 154         const BASE Aij_real = CONST_REAL(A, lda * j + i);
 155         const BASE Aij_imag = conj * CONST_IMAG(A, lda * j + i);
 156         const BASE x_real = REAL(X, jx);
 157         const BASE x_imag = IMAG(X, jx);
 158         tmp_real -= Aij_real * x_real - Aij_imag * x_imag;
 159         tmp_imag -= Aij_real * x_imag + Aij_imag * x_real;
 160         jx += incX;
 161       }
 162       if (nonunit) {
 163         const BASE a_real = CONST_REAL(A, lda * i + i);
 164         const BASE a_imag = conj * CONST_IMAG(A, lda * i + 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 if ((order == CblasRowMajor && Trans == CblasTrans && Uplo == CblasLower)
 177              || (order == CblasColMajor && Trans == CblasNoTrans && Uplo == CblasUpper)) {
 178
 179     /* backsubstitution */
 180
 181     ix = OFFSET(N, incX) + incX * (N - 1);
 182
 183     if (nonunit) {
 184       const BASE a_real = CONST_REAL(A, lda * (N - 1) + (N - 1));
 185       const BASE a_imag = conj * CONST_IMAG(A, lda * (N - 1) + (N - 1));
 186       const BASE x_real = REAL(X, ix);
 187       const BASE x_imag = IMAG(X, ix);
 188       const BASE s = xhypot(a_real, a_imag);
 189       const BASE b_real = a_real / s;
 190       const BASE b_imag = a_imag / s;
 191       REAL(X, ix) = (x_real * b_real + x_imag * b_imag) / s;
 192       IMAG(X, ix) = (x_imag * b_real - b_imag * x_real) / s;
 193     }
 194
 195     ix -= incX;
 196
 197     for (i = N - 1; i > 0 && i--;) {
 198       BASE tmp_real = REAL(X, ix);
 199       BASE tmp_imag = IMAG(X, ix);
 200       jx = ix + incX;
 201       for (j = i + 1; j < N; j++) {
 202         const BASE Aij_real = CONST_REAL(A, lda * j + i);
 203         const BASE Aij_imag = conj * CONST_IMAG(A, lda * j + i);
 204         const BASE x_real = REAL(X, jx);
 205         const BASE x_imag = IMAG(X, jx);
 206         tmp_real -= Aij_real * x_real - Aij_imag * x_imag;
 207         tmp_imag -= Aij_real * x_imag + Aij_imag * x_real;
 208         jx += incX;
 209       }
 210
 211       if (nonunit) {
 212         const BASE a_real = CONST_REAL(A, lda * i + i);
 213         const BASE a_imag = conj * CONST_IMAG(A, lda * i + i);
 214         const BASE s = xhypot(a_real, a_imag);
 215         const BASE b_real = a_real / s;
 216         const BASE b_imag = a_imag / s;
 217         REAL(X, ix) = (tmp_real * b_real + tmp_imag * b_imag) / s;
 218         IMAG(X, ix) = (tmp_imag * b_real - tmp_real * b_imag) / s;
 219       } else {
 220         REAL(X, ix) = tmp_real;
 221         IMAG(X, ix) = tmp_imag;
 222       }
 223       ix -= incX;
 224     }
 225   } else {
 226     BLAS_ERROR("unrecognized operation");
 227   }
 228 }