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 } 180