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