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

   1 /* blas/source_spmv.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
  23
  24   if (alpha == 0.0 && beta == 1.0)
  25     return;
  26
  27   /* form  y := beta*y */
  28   if (beta == 0.0) {
  29     INDEX iy = OFFSET(N, incY);
  30     for (i = 0; i < N; i++) {
  31       Y[iy] = 0.0;
  32       iy += incY;
  33     }
  34   } else if (beta != 1.0) {
  35     INDEX iy = OFFSET(N, incY);
  36     for (i = 0; i < N; i++) {
  37       Y[iy] *= beta;
  38       iy += incY;
  39     }
  40   }
  41
  42   if (alpha == 0.0)
  43     return;
  44
  45   /* form  y := alpha*A*x + y */
  46
  47   if ((order == CblasRowMajor && Uplo == CblasUpper)
  48       || (order == CblasColMajor && Uplo == CblasLower)) {
  49     INDEX ix = OFFSET(N, incX);
  50     INDEX iy = OFFSET(N, incY);
  51     for (i = 0; i < N; i++) {
  52       BASE tmp1 = alpha * X[ix];
  53       BASE tmp2 = 0.0;
  54       const INDEX j_min = i + 1;
  55       const INDEX j_max = N;
  56       INDEX jx = OFFSET(N, incX) + j_min * incX;
  57       INDEX jy = OFFSET(N, incY) + j_min * incY;
  58
  59       Y[iy] += tmp1 * Ap[TPUP(N, i, i)];
  60
  61       for (j = j_min; j < j_max; j++) {
  62         const BASE apk = Ap[TPUP(N, i, j)];
  63         Y[jy] += tmp1 * apk;
  64         tmp2 += apk * X[jx];
  65         jy += incY;
  66         jx += incX;
  67       }
  68       Y[iy] += alpha * tmp2;
  69       ix += incX;
  70       iy += incY;
  71     }
  72   } else if ((order == CblasRowMajor && Uplo == CblasLower)
  73              || (order == CblasColMajor && Uplo == CblasUpper)) {
  74     INDEX ix = OFFSET(N, incX);
  75     INDEX iy = OFFSET(N, incY);
  76     for (i = 0; i < N; i++) {
  77       BASE tmp1 = alpha * X[ix];
  78       BASE tmp2 = 0.0;
  79
  80       const INDEX j_min = 0;
  81       const INDEX j_max = i;
  82       INDEX jx = OFFSET(N, incX) + j_min * incX;
  83       INDEX jy = OFFSET(N, incY) + j_min * incY;
  84
  85       Y[iy] += tmp1 * Ap[TPLO(N, i, i)];
  86
  87       for (j = j_min; j < j_max; j++) {
  88         const BASE apk = Ap[TPLO(N, i, j)];
  89         Y[jy] += tmp1 * apk;
  90         tmp2 += apk * X[jx];
  91         jy += incY;
  92         jx += incX;
  93       }
  94       Y[iy] += alpha * tmp2;
  95       ix += incX;
  96       iy += incY;
  97     }
  98   } else {
  99     BLAS_ERROR("unrecognized operation");
 100   }
 101 }