ifp_spmm.cpp

00001 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00002 /*             ********   ***                                 SparseLib++    */
00003 /*          *******  **  ***       ***      ***               v. 1.5         */
00004 /*           *****      ***     ******** ********                            */
00005 /*            *****    ***     ******** ********              R. Pozo        */
00006 /*       **  *******  ***   **   ***      ***                 K. Remington   */
00007 /*        ********   ********                                 A. Lumsdaine   */
00008 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00009 /*                                                                           */
00010 /*                                                                           */
00011 /*                     SparseLib++ : Sparse Matrix Library                   */
00012 /*                                                                           */
00013 /*               National Institute of Standards and Technology              */
00014 /*                        University of Notre Dame                           */
00015 /*              Authors: R. Pozo, K. Remington, A. Lumsdaine                 */
00016 /*                                                                           */
00017 /*                                 NOTICE                                    */
00018 /*                                                                           */
00019 /* Permission to use, copy, modify, and distribute this software and         */
00020 /* its documentation for any purpose and without fee is hereby granted       */
00021 /* provided that the above notice appear in all copies and supporting        */
00022 /* documentation.                                                            */
00023 /*                                                                           */
00024 /* Neither the Institutions (National Institute of Standards and Technology, */
00025 /* University of Notre Dame) nor the Authors make any representations about  */
00026 /* the suitability of this software for any purpose.  This software is       */
00027 /* provided ``as is'' without expressed or implied warranty.                 */
00028 /*                                                                           */
00029 /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
00030 
00031 /*
00032  * Home Grown Sparse BLAS
00033  *
00034  * These are just a subset of the functions described in SPARKER
00035  * Working Note #3.
00036  *
00037  * Would be great if these could be templated some day
00038  *
00039  */
00040 
00041 #include <stdlib.h>
00042 #include <iostream>
00043 using namespace std;
00044 #include "ifp_spblas.h"
00045 
00046 #define _SpMatVal(_a,_lda,_row,_col) ((_a)[(_lda)*(_col)+(_row)])
00047 
00048 static void CoordMatVec_float(int m, int n, int k, const float &alpha,
00049         const float *val, const int *indx, const int *jndx,
00050         const int &nnz,
00051         const float *b, int ldb, float *c, int ldc)
00052 {
00053   int i, j;
00054 
00055   // To make the compiler happy
00056   if (k && m)
00057     ;
00058 
00059   // Frob these so we can use one-based indexing externally
00060   b -= 1;
00061   c -= 1;
00062 
00063   if (alpha == 1.0) {
00064     if (n == 1)
00065       for (j = 0; j < nnz; j++)
00066     c[indx[j]] += b[jndx[j]] * val[j];
00067     else
00068       for (i = 0; i < n; i++)
00069     for (j = 0; j < nnz; j++)
00070       _SpMatVal(c, ldc, indx[j], i) += _SpMatVal(b, ldb, indx[j], i) * val[j];
00071   } else {
00072     if (n == 1)
00073       for (j = 0; j < nnz; j++)
00074     c[indx[j]] += alpha * b[jndx[j]] * val[j];
00075     else
00076       for (i = 0; i < n; i++)
00077     for (j = 0; j < nnz; j++)
00078       _SpMatVal(c, ldc, indx[j], i) +=
00079         alpha * _SpMatVal(b, ldb, indx[j], i) * val[j];
00080   }
00081 }
00082 
00083 static void CoordMatVec_double(int m, int n, int k, const double &alpha,
00084         const double *val, const int *indx, const int *jndx,
00085         const int &nnz,
00086         const double *b, int ldb, double *c, int ldc)
00087 {
00088   int i, j;
00089 
00090   // To make the compiler happy
00091   if (k && m)
00092     ;
00093 
00094   // Frob these so we can use one-based indexing externally
00095   b -= 1;
00096   c -= 1;
00097 
00098   if (alpha == 1.0) {
00099     if (n == 1)
00100       for (j = 0; j < nnz; j++)
00101     c[indx[j]] += b[jndx[j]] * val[j];
00102     else
00103       for (i = 0; i < n; i++)
00104     for (j = 0; j < nnz; j++)
00105       _SpMatVal(c, ldc, indx[j], i) += _SpMatVal(b, ldb, indx[j], i) * val[j];
00106   } else {
00107     if (n == 1)
00108       for (j = 0; j < nnz; j++)
00109     c[indx[j]] += alpha * b[jndx[j]] * val[j];
00110     else
00111       for (i = 0; i < n; i++)
00112     for (j = 0; j < nnz; j++)
00113       _SpMatVal(c, ldc, indx[j], i) +=
00114         alpha * _SpMatVal(b, ldb, indx[j], i) * val[j];
00115   }
00116 }
00117 
00118 static void
00119 CompColMatVec_double(int m, int n, int k, const double &alpha,
00120             const double *val, const int *indx, const int *pntr,
00121             const double *b, int ldb, double *c, int ldc)
00122 {
00123   int i, j, l;
00124 
00125   if (alpha == 0.0)
00126     return;
00127 
00128   // To make the compiler happy
00129   if (m)
00130     ;
00131 
00132   // Frob these so we can use one-based indexing externally
00133   c -= 1;
00134   val -= pntr[0];
00135   indx -= pntr[0];
00136 
00137   if (alpha == 1.0) {
00138     if (n == 1)
00139       for (i = 0; i < k; i++)
00140     for (j = pntr[i]; j < pntr[i+1]; j++)
00141       c[indx[j]] += b[i] * val[j];
00142     else
00143       for (l = 0; l < n; l++)
00144     for (i = 0; i < k; i++)
00145       for (j = pntr[i]; j < pntr[i+1]; j++)
00146         _SpMatVal(c, ldc, indx[j], l) += _SpMatVal(b, ldb, i, l) * val[j];
00147   } else {
00148     if (n == 1)
00149       for (i = 0; i < k; i++)
00150     for (j = pntr[i]; j < pntr[i+1]; j++)
00151       c[indx[j]] += alpha * b[i] * val[j];
00152     else
00153       for (l = 0; l < n; l++)
00154     for (i = 0; i < k; i++)
00155       for (j = pntr[i]; j < pntr[i+1]; j++)
00156         _SpMatVal(c, ldc, indx[j], l) +=
00157           alpha * _SpMatVal(b, ldb, i, l) * val[j];
00158   }
00159 }
00160 
00161 static void CompColMatVec_float(int m, int n, int k, const float &alpha,
00162             const float *val, const int *indx, const int *pntr,
00163             const float *b, int ldb, float *c, int ldc)
00164 {
00165   int i, j, l;
00166 
00167   if (alpha == 0.0)
00168     return;
00169 
00170   // To make the compiler happy
00171   if (m)
00172     ;
00173 
00174   // Frob these so we can use one-based indexing externally
00175   c -= 1;
00176   val -= pntr[0];
00177   indx -= pntr[0];
00178 
00179   if (alpha == 1.0) {
00180     if (n == 1)
00181       for (i = 0; i < k; i++)
00182     for (j = pntr[i]; j < pntr[i+1]; j++)
00183       c[indx[j]] += b[i] * val[j];
00184     else
00185       for (l = 0; l < n; l++)
00186     for (i = 0; i < k; i++)
00187       for (j = pntr[i]; j < pntr[i+1]; j++)
00188         _SpMatVal(c, ldc, indx[j], l) += _SpMatVal(b, ldb, i, l) * val[j];
00189   } else {
00190     if (n == 1)
00191       for (i = 0; i < k; i++)
00192     for (j = pntr[i]; j < pntr[i+1]; j++)
00193       c[indx[j]] += alpha * b[i] * val[j];
00194     else
00195       for (l = 0; l < n; l++)
00196     for (i = 0; i < k; i++)
00197       for (j = pntr[i]; j < pntr[i+1]; j++)
00198         _SpMatVal(c, ldc, indx[j], l) +=
00199           alpha * _SpMatVal(b, ldb, i, l) * val[j];
00200   }
00201 }
00202 
00203 
00204 static void
00205 CompRowMatVec_double(int m, int n, int k, const double &alpha,
00206             const double *val, const int *indx, const int *pntr,
00207             const double *b, int ldb, double *c, int ldc)
00208 {
00209   int i, j, l;
00210 
00211   if (alpha == 0.0)
00212     return;
00213 
00214   // To make the compiler happy
00215   if (m || k)
00216     ;
00217 
00218   // Frob these so we can use one-based indexing externally
00219   b -= 1;
00220   val -= pntr[0];
00221   indx -= pntr[0];
00222 
00223   if (alpha == 1.0) {
00224     if (n == 1)
00225       for (i = 0; i < m; i++)
00226     for (j = pntr[i]; j < pntr[i+1]; j++)
00227       c[i] += b[indx[j]] * val[j];
00228     else
00229       for (l = 0; l < n; l++)
00230     for (i = 0; i < m; i++)
00231       for (j = pntr[i]; j < pntr[i+1]; j++)
00232         _SpMatVal(c, ldc, i, l) += _SpMatVal(b, ldb, indx[j], l) * val[j];
00233   } else {
00234     if (n == 1)
00235       for (i = 0; i < m; i++)
00236     for (j = pntr[i]; j < pntr[i+1]; j++)
00237       c[i] += alpha * b[indx[j]] * val[j];
00238     else
00239       for (l = 0; l < n; l++)
00240     for (i = 0; i < m; i++)
00241       for (j = pntr[i]; j < pntr[i+1]; j++)
00242         _SpMatVal(c, ldc, i, l) +=
00243           alpha * _SpMatVal(b, ldb, indx[j], l) * val[j];
00244   }
00245 }
00246 
00247 
00248 static void
00249 CompRowMatVec_float(int m, int n, int k, const float &alpha,
00250             const float *val, const int *indx, const int *pntr,
00251             const float *b, int ldb, float *c, int ldc)
00252 {
00253   int i, j, l;
00254 
00255   if (alpha == 0.0)
00256     return;
00257 
00258   // To make the compiler happy
00259   if (m || k)
00260     ;
00261 
00262   // Frob these so we can use one-based indexing externally
00263   b -= 1;
00264   val -= pntr[0];
00265   indx -= pntr[0];
00266 
00267   if (alpha == 1.0) {
00268     if (n == 1)
00269       for (i = 0; i < m; i++)
00270     for (j = pntr[i]; j < pntr[i+1]; j++)
00271       c[i] += b[indx[j]] * val[j];
00272     else
00273       for (l = 0; l < n; l++)
00274     for (i = 0; i < m; i++)
00275       for (j = pntr[i]; j < pntr[i+1]; j++)
00276         _SpMatVal(c, ldc, i, l) += _SpMatVal(b, ldb, indx[j], l) * val[j];
00277   } else {
00278     if (n == 1)
00279       for (i = 0; i < m; i++)
00280     for (j = pntr[i]; j < pntr[i+1]; j++)
00281       c[i] += alpha * b[indx[j]] * val[j];
00282     else
00283       for (l = 0; l < n; l++)
00284     for (i = 0; i < m; i++)
00285       for (j = pntr[i]; j < pntr[i+1]; j++)
00286         _SpMatVal(c, ldc, i, l) +=
00287           alpha * _SpMatVal(b, ldb, indx[j], l) * val[j];
00288   }
00289 }
00290 
00291 static void
00292 ScaleRectangularArray_double(int m, int n, double *c, int ldc, 
00293     const double &beta)
00294 {
00295   int i, j;
00296 
00297   if (beta == 1.0)
00298     return;
00299 
00300   if (beta == 0.0) {
00301     if (n == 1)
00302       for (j = 0; j < m; j++)
00303     c[j] = 0.0;
00304     else
00305       for (i = 0; i < n; i++)
00306     for (j = 0; j < m; j++)
00307       _SpMatVal(c, ldc, j, i) = 0.0;
00308   } else {
00309     if (n == 1)
00310       for (j = 0; j < m; j++)
00311     c[j] *= beta;
00312     else
00313       for (i = 0; i < n; i++)
00314     for (j = 0; j < m; j++)
00315       _SpMatVal(c, ldc, j, i) *= beta;
00316   }
00317 }
00318 
00319 static void
00320 ScaleRectangularArray_float(int m, int n, float *c, int ldc, 
00321     const double &beta)
00322 {
00323   int i, j;
00324 
00325   if (beta == 1.0)
00326     return;
00327 
00328   if (beta == 0.0) {
00329     if (n == 1)
00330       for (j = 0; j < m; j++)
00331     c[j] = 0.0;
00332     else
00333       for (i = 0; i < n; i++)
00334     for (j = 0; j < m; j++)
00335       _SpMatVal(c, ldc, j, i) = 0.0;
00336   } else {
00337     if (n == 1)
00338       for (j = 0; j < m; j++)
00339     c[j] *= beta;
00340     else
00341       for (i = 0; i < n; i++)
00342     for (j = 0; j < m; j++)
00343       _SpMatVal(c, ldc, j, i) *= beta;
00344   }
00345 }
00346 
00347 /*
00348  * dcoom -- coordinate format matrix-matrix multiply
00349  *
00350  * C <- alpha A B + beta C
00351  *
00352  * Arguments:
00353  *
00354  * int &transa  Indicates how to operate with the sparse matrix
00355  *      0 : operate with matrix
00356  *      1 : operate with transpose matrix
00357  *      2 : operate with conjugate transpose matrix
00358  *
00359  * int &m   Number of rows in matrix c
00360  *
00361  * int &n   Number of columns in matrix c
00362  *
00363  * int &k   Number of rows in matrix b
00364  *
00365  * double &alpha Scalar parameter
00366  *
00367  * double &beta Scalar parameter
00368  *
00369  * int descra[] Descriptor argument.  Nine element integer array
00370  *      descra[0] matrix structure
00371  *          0 : general
00372  *          1 : symmetric
00373  *          2 : Hermition
00374  *          3 : Triangular
00375  *          4 : Anti-Symmetric
00376  *          5 : Diagonal
00377  *      descra[1] upper/lower triangular indicator
00378  *          1 : lower
00379  *          2 : upper
00380  *      descra[2] main diagonal type
00381  *          0 : non-unit
00382  *          1 : unit
00383  *      descra[4] repeated indices?
00384  *          0 : unknown
00385  *          1 : no repeated indices
00386  *
00387  *
00388  * double *val  scalar array of length nnz containing matrix entries
00389  *
00390  * int *indx    integer array of length nnz containing row indices
00391  *
00392  * int *jndx    integer array of length nnz containing column indices
00393  *
00394  * double *b    rectangular array with first dimension ldb
00395  *
00396  * double *c    rectangular array with first dimension ldc
00397  *
00398  * double *work scratch array of length lwork.  lwork should be at least
00399  *      max(m,n)
00400  *
00401  */
00402 void F77NAME(scoomm)
00403   (const int &transa, const int &m, const int &n, const int &k,
00404    const float &alpha,
00405    const int descra[], const float *val,
00406    const int *indx, const int *jndx, const int &nnz,
00407    const float *b, const int &ldb,
00408    const float &beta, float *c, const int &ldc,
00409    float *work, const int &lwork)
00410 {
00411   if (descra[0] != 0) {
00412     cerr << "Must have general matrix" << endl;
00413     exit(1);
00414   }
00415 
00416   // To make the compiler happy
00417   if (work && lwork)
00418     ;
00419 
00420   ScaleRectangularArray_float(m, n, c, ldc, beta);
00421 
00422   if (alpha == 0.0)
00423     return;
00424 
00425   // Use this hack if transpose is desired
00426   if (transa == 1 || transa == 2) {
00427     const int *itmp = indx;
00428     indx = jndx;
00429     jndx = itmp;
00430   }
00431   CoordMatVec_float(m, n, k, alpha, val, indx, jndx, nnz, b, ldb, c, ldc);
00432 }
00433 
00434 
00435 void F77NAME(dcoomm)
00436   (const int &transa, const int &m, const int &n, const int &k,
00437    const double &alpha,
00438    const int descra[], const double *val,
00439    const int *indx, const int *jndx, const int &nnz,
00440    const double *b, const int &ldb,
00441    const double &beta, double *c, const int &ldc,
00442    double *work, const int &lwork)
00443 {
00444   if (descra[0] != 0) {
00445     cerr << "Must have general matrix" << endl;
00446     exit(1);
00447   }
00448 
00449   // To make the compiler happy
00450   if (work && lwork)
00451     ;
00452 
00453   ScaleRectangularArray_double(m, n, c, ldc, beta);
00454 
00455   if (alpha == 0.0)
00456     return;
00457 
00458   // Use this hack if transpose is desired
00459   if (transa == 1 || transa == 2) {
00460     const int *itmp = indx;
00461     indx = jndx;
00462     jndx = itmp;
00463   }
00464   CoordMatVec_double(m, n, k, alpha, val, indx, jndx, nnz, b, ldb, c, ldc);
00465 }
00466 
00467 
00468 
00469 
00470 /*
00471  * dcscm -- comp sparse column matrix-matrix multiply
00472  *
00473  * Arguments:
00474  *
00475  * int &transa  Indicates how to operate with the sparse matrix
00476  *      0 : operate with matrix
00477  *      1 : operate with transpose matrix
00478  *      2 : operate with conjugate transpose matrix
00479  *
00480  * int &m   Number of rows in matrix c
00481  *
00482  * int &n   Number of columns in matrix c
00483  *
00484  * int &k   Number of rows in matrix b
00485  *
00486  * double &alpha Scalar parameter
00487  *
00488  * double &beta Scalar parameter
00489  *
00490  * int descra[] Descriptor argument.  Nine element integer array
00491  *      descra[0] matrix structure
00492  *          0 : general
00493  *          1 : symmetric
00494  *          2 : Hermition
00495  *          3 : Triangular
00496  *          4 : Anti-Symmetric
00497  *          5 : Diagonal
00498  *      descra[1] upper/lower triangular indicator
00499  *          1 : lower
00500  *          2 : upper
00501  *      descra[2] main diagonal type
00502  *          0 : non-unit
00503  *          1 : unit
00504  *
00505  * double *val  scalar array of length nnz containing matrix entries
00506  *
00507  * int *indx    integer array of length nnz containing row indices
00508  *
00509  * int *pntr    integer array of length k+1 such that pntr(j)-pntr(1)
00510  *      points to location in val of the first element in column j
00511  *
00512  * double *b    rectangular array with first dimension ldb
00513  *
00514  * double *c    rectangular array with first dimension ldc
00515  *
00516  * double *work scratch array of length lwork.  lwork should be at least
00517  *      max(m,n)
00518  *
00519  */
00520 void F77NAME(scscmm)
00521   (const int &transa, const int &m, const int &n, const int &k,
00522    const float &alpha,
00523    const int descra[], const float *val,
00524    const int *indx, const int *pntr, const float *b, int &ldb,
00525    const float &beta, float *c, const int &ldc,
00526    float *work, const int &lwork)
00527 {
00528   if (descra[0] != 0) {
00529     cerr << "Must have general matrix" << endl;
00530     exit(1);
00531   }
00532 
00533   // To make the compiler happy
00534   if (work && lwork)
00535     ;
00536 
00537   ScaleRectangularArray_float(m, n, c, ldc, beta);
00538 
00539   if (transa == 1 || transa == 2)
00540     CompRowMatVec_float(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
00541   else
00542     CompColMatVec_float(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
00543 }
00544 
00545 
00546 void F77NAME(dcscmm)
00547   (const int &transa, const int &m, const int &n, const int &k,
00548    const double &alpha,
00549    const int descra[], const double *val,
00550    const int *indx, const int *pntr, const double *b, int &ldb,
00551    const double &beta, double *c, const int &ldc,
00552    double *work, const int &lwork)
00553 {
00554   if (descra[0] != 0) {
00555     cerr << "Must have general matrix" << endl;
00556     exit(1);
00557   }
00558 
00559   // To make the compiler happy
00560   if (work && lwork)
00561     ;
00562 
00563   ScaleRectangularArray_double(m, n, c, ldc, beta);
00564 
00565   if (transa == 1 || transa == 2)
00566     CompRowMatVec_double(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
00567   else
00568     CompColMatVec_double(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
00569 }
00570 
00571 
00572 
00573 
00574 /*
00575  * dcsrm -- comp sparse row matrix-matrix multiply
00576  *
00577  * Arguments:
00578  *
00579  * int &transa  Indicates how to operate with the sparse matrix
00580  *      0 : operate with matrix
00581  *      1 : operate with transpose matrix
00582  *      2 : operate with conjugate transpose matrix
00583  *
00584  * int &m   Number of rows in matrix c
00585  *
00586  * int &n   Number of columns in matrix c
00587  *
00588  * int &k   Number of rows in matrix b
00589  *
00590  * double &alpha Scalar parameter
00591  *
00592  * double &beta Scalar parameter
00593  *
00594  * int descra[] Descriptor argument.  Nine element integer array
00595  *      descra[0] matrix structure
00596  *          0 : general
00597  *          1 : symmetric
00598  *          2 : Hermition
00599  *          3 : Triangular
00600  *          4 : Anti-Symmetric
00601  *          5 : Diagonal
00602  *      descra[1] upper/lower triangular indicator
00603  *          1 : lower
00604  *          2 : upper
00605  *      descra[2] main diagonal type
00606  *          0 : non-unit
00607  *          1 : unit
00608  *
00609  * double *val  scalar array of length nnz containing matrix entries
00610  *
00611  * int *indx    integer array of length nnz containing column indices
00612  *
00613  * int *pntr    integer array of length k+1 such that pntr(j)-pntr(1)
00614  *      points to location in val of the first element in row j
00615  *
00616  * double *b    rectangular array with first dimension ldb
00617  *
00618  * double *c    rectangular array with first dimension ldc
00619  *
00620  * double *work scratch array of length lwork.  lwork should be at least
00621  *      max(m,n)
00622  *
00623  */
00624 void F77NAME(scsrmm)
00625   (const int &transa, const int &m, const int &n, const int &k,
00626    const float &alpha,
00627    const int descra[], const float *val,
00628    const int *indx, const int *pntr, const float *b, int &ldb,
00629    const float &beta, float *c, const int &ldc,
00630    float *work, const int &lwork)
00631 {
00632   if (descra[0] != 0) {
00633     cerr << "Must have general matrix" << endl;
00634     exit(1);
00635   }
00636 
00637   // To make the compiler happy
00638   if (work && lwork)
00639     ;
00640 
00641   ScaleRectangularArray_float(m, n, c, ldc, beta);
00642 
00643   if (transa == 1 || transa == 2)
00644     CompColMatVec_float(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
00645   else
00646     CompRowMatVec_float(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
00647 }
00648 
00649 
00650 void F77NAME(dcsrmm)
00651   (const int &transa, const int &m, const int &n, const int &k,
00652    const double &alpha,
00653    const int descra[], const double *val,
00654    const int *indx, const int *pntr, const double *b, int &ldb,
00655    const double &beta, double *c, const int &ldc,
00656    double *work, const int &lwork)
00657 {
00658   if (descra[0] != 0) {
00659     cerr << "Must have general matrix" << endl;
00660     exit(1);
00661   }
00662 
00663   // To make the compiler happy
00664   if (work && lwork)
00665     ;
00666 
00667   ScaleRectangularArray_double(m, n, c, ldc, beta);
00668 
00669   if (transa == 1 || transa == 2)
00670     CompColMatVec_double(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
00671   else
00672     CompRowMatVec_double(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
00673 }
00674 
00675 
00676 
00677 
00678 
 All Classes Files Functions Variables Enumerations Friends
Generated on Wed Apr 13 10:05:22 2011 for IFPACK by  doxygen 1.6.3