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

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