1 /* blas/source_trmv_r.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 INDEX i, j; 22 const int nonunit = (Diag == CblasNonUnit); 23 const int Trans = (TransA != CblasConjTrans) ? TransA : CblasTrans; 24 25 if ((order == CblasRowMajor && Trans == CblasNoTrans && Uplo == CblasUpper) 26 || (order == CblasColMajor && Trans == CblasTrans && Uplo == CblasLower)) { 27 /* form x := A*x */ 28 29 INDEX ix = OFFSET(N, incX); 30 for (i = 0; i < N; i++) { 31 BASE temp = 0.0; 32 const INDEX j_min = i + 1; 33 const INDEX j_max = N; 34 INDEX jx = OFFSET(N, incX) + j_min * incX; 35 for (j = j_min; j < j_max; j++) { 36 temp += X[jx] * A[lda * i + j]; 37 jx += incX; 38 } 39 if (nonunit) { 40 X[ix] = temp + X[ix] * A[lda * i + i]; 41 } else { 42 X[ix] += temp; 43 } 44 ix += incX; 45 } 46 } else if ((order == CblasRowMajor && Trans == CblasNoTrans && Uplo == CblasLower) 47 || (order == CblasColMajor && Trans == CblasTrans && Uplo == CblasUpper)) { 48 INDEX ix = OFFSET(N, incX) + (N - 1) * incX; 49 for (i = N; i > 0 && i--;) { 50 BASE temp = 0.0; 51 const INDEX j_min = 0; 52 const INDEX j_max = i; 53 INDEX jx = OFFSET(N, incX) + j_min * incX; 54 for (j = j_min; j < j_max; j++) { 55 temp += X[jx] * A[lda * i + j]; 56 jx += incX; 57 } 58 if (nonunit) { 59 X[ix] = temp + X[ix] * A[lda * i + i]; 60 } else { 61 X[ix] += temp; 62 } 63 ix -= incX; 64 } 65 } else if ((order == CblasRowMajor && Trans == CblasTrans && Uplo == CblasUpper) 66 || (order == CblasColMajor && Trans == CblasNoTrans && Uplo == CblasLower)) { 67 /* form x := A'*x */ 68 INDEX ix = OFFSET(N, incX) + (N - 1) * incX; 69 for (i = N; i > 0 && i--;) { 70 BASE temp = 0.0; 71 const INDEX j_min = 0; 72 const INDEX j_max = i; 73 INDEX jx = OFFSET(N, incX) + j_min * incX; 74 for (j = j_min; j < j_max; j++) { 75 temp += X[jx] * A[lda * j + i]; 76 jx += incX; 77 } foldSelectWithBinaryOp(Value * Cond,Value * TrueVal,Value * FalseVal)78 if (nonunit) { 79 X[ix] = temp + X[ix] * A[lda * i + i]; 80 } else { 81 X[ix] += temp; 82 } 83 ix -= incX; 84 } 85 } else if ((order == CblasRowMajor && Trans == CblasTrans && Uplo == CblasLower) 86 || (order == CblasColMajor && Trans == CblasNoTrans && Uplo == CblasUpper)) { 87 INDEX ix = OFFSET(N, incX); 88 for (i = 0; i < N; i++) { 89 BASE temp = 0.0; 90 const INDEX j_min = i + 1; 91 const INDEX j_max = N; 92 INDEX jx = OFFSET(N, incX) + (i + 1) * incX; 93 for (j = j_min; j < j_max; j++) { 94 temp += X[jx] * A[lda * j + i]; 95 jx += incX; 96 } 97 if (nonunit) { 98 X[ix] = temp + X[ix] * A[lda * i + i]; 99 } else { 100 X[ix] += temp; 101 } 102 ix += incX; 103 } 104 } else { 105 BLAS_ERROR("unrecognized operation"); 106 } 107 } 108