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

   1 /* blas/source_tpmv_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 (N == 0)
  26     return;
  27
  28   if ((order == CblasRowMajor && Trans == CblasNoTrans && Uplo == CblasUpper)
  29       || (order == CblasColMajor && Trans == CblasTrans && Uplo == CblasLower)) {
  30     /* form  x:= A*x */
  31
  32     INDEX ix = OFFSET(N, incX);
  33     for (i = 0; i < N; i++) {
  34       BASE atmp = Ap[TPUP(N, i, i)];
  35       BASE temp = (nonunit ? X[ix] * atmp : X[ix]);
  36       INDEX jx = OFFSET(N, incX) + (i + 1) * incX;
  37       for (j = i + 1; j < N; j++) {
  38         atmp = Ap[TPUP(N, i, j)];
  39         temp += atmp * X[jx];
  40         jx += incX;
  41       }
  42       X[ix] = temp;
  43       ix += incX;
  44     }
  45   } else if ((order == CblasRowMajor && Trans == CblasNoTrans && Uplo == CblasLower)
  46              || (order == CblasColMajor && Trans == CblasTrans && Uplo == CblasUpper)) {
  47
  48     INDEX ix = OFFSET(N, incX) + (N - 1) * incX;
  49     for (i = N; i > 0 && i--;) {
  50       BASE atmp = Ap[TPLO(N, i, i)];
  51       BASE temp = (nonunit ? X[ix] * atmp : X[ix]);
  52       INDEX jx = OFFSET(N, incX);
  53       for (j = 0; j < i; j++) {
  54         atmp = Ap[TPLO(N, i, j)];
  55         temp += atmp * X[jx];
  56         jx += incX;
  57       }
  58       X[ix] = temp;
  59       ix -= incX;
  60     }
  61   } else if ((order == CblasRowMajor && Trans == CblasTrans && Uplo == CblasUpper)
  62              || (order == CblasColMajor && Trans == CblasNoTrans && Uplo == CblasLower)) {
  63     /* form  x := A'*x */
  64
  65     INDEX ix = OFFSET(N, incX) + (N - 1) * incX;
  66     for (i = N; i > 0 && i--;) {
  67       BASE atmp = Ap[TPUP(N, i, i)];
  68       BASE temp = (nonunit ? X[ix] * atmp : X[ix]);
  69       INDEX jx = OFFSET(N, incX);
  70       for (j = 0; j < i; j++) {
  71         atmp = Ap[TPUP(N, j, i)];
  72         temp += atmp * X[jx];
  73         jx += incX;
  74       }
  75       X[ix] = temp;
  76       ix -= incX;
  77     }
  78   } else if ((order == CblasRowMajor && Trans == CblasTrans && Uplo == CblasLower)
  79              || (order == CblasColMajor && Trans == CblasNoTrans && Uplo == CblasUpper)) {
  80
  81     INDEX ix = OFFSET(N, incX);
  82     for (i = 0; i < N; i++) {
  83       BASE atmp = Ap[TPLO(N, i, i)];
  84       BASE temp = (nonunit ? X[ix] * atmp : X[ix]);
  85       INDEX jx = OFFSET(N, incX) + (i + 1) * incX;
  86       for (j = i + 1; j < N; j++) {
  87         atmp = Ap[TPLO(N, j, i)];
  88         temp += atmp * X[jx];
  89         jx += incX;
  90       }
  91       X[ix] = temp;
  92       ix += incX;
  93     }
  94   } else {
  95     BLAS_ERROR("unrecognized operation");
  96   }
  97 }