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