Ifpack Package Browser (Single Doxygen Collection) Development
ifp_spsm.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 
00051 /*
00052  * int &m   Number of rows in matrix c
00053  *
00054  * int &n   Number of columns in matrix c
00055  *
00056  * int &k   Number of rows in matrix b
00057  *
00058  * Assume diagonal elements are in proper place
00059  *
00060  * unitd = 1    D = I
00061  * unitd = 2    left (row scaling)
00062  * unitd = 3    right (column scaling)
00063  */
00064 
00065 static void
00066 CopyRectangularArray_double(int m, int n,
00067              const double *b, int ldb, double *c, int ldc)
00068 {
00069   int i, l;
00070 
00071   if (b == c)
00072     return;
00073 
00074   if (n == 1)
00075     for (i = 0; i < m; i++)
00076       c[i] = b[i];
00077   else
00078     for (l = 0; l < n; l++)
00079       for (i = 0; i < m; i++)
00080     _SpMatVal(c, ldc, i, l) = _SpMatVal(b, ldb, i, l);
00081 }
00082 
00083 
00084 static void
00085 CopyRectangularArray_float(int m, int n,
00086              const float *b, int ldb, float *c, int ldc)
00087 {
00088   int i, l;
00089 
00090   if (b == c)
00091     return;
00092 
00093   if (n == 1)
00094     for (i = 0; i < m; i++)
00095       c[i] = b[i];
00096   else
00097     for (l = 0; l < n; l++)
00098       for (i = 0; i < m; i++)
00099     _SpMatVal(c, ldc, i, l) = _SpMatVal(b, ldb, i, l);
00100 }
00101 
00102 static void
00103 CompCol_LowerUnitSolve_double(int m, int n, int unitd, const double *dv, 
00104       double alpha, const double *val, const int *indx, const int *pntr,
00105       const double *b, int ldb, double *c, int ldc)
00106 {
00107   int i, j, l;
00108   double z;
00109 
00110   // To make the compiler happy
00111   if (dv)
00112     ;
00113 
00114   if (unitd != 1) {
00115     cerr << "unitd != 1 not implemented" << endl;
00116     exit(1);
00117   }
00118 
00119   if (alpha == 0.0)
00120     return;
00121 
00122   CopyRectangularArray_double(m, n, b, ldb, &c[pntr[0]-1], ldc);
00123 
00124   c -= 1;
00125   val -= pntr[0];
00126   indx -= pntr[0];
00127 
00128   if (alpha == 1.0) {
00129     if (n == 1)
00130       for (i = 0; i < m; i++) {
00131     z = c[i+pntr[0]];
00132     for (j = pntr[i]; j < pntr[i+1]; j++)
00133       c[indx[j]] -= z * val[j];
00134       }
00135     else
00136       for (l = 0; l < n; l++)
00137     for (i = 0; i < m; i++) {
00138       z = _SpMatVal(c, ldc, i+pntr[0], l);
00139       for (j = pntr[i]; j < pntr[i+1]; j++)
00140         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00141     }
00142   } else {
00143     if (n == 1)
00144       for (i = 0; i < m; i++) {
00145     z = alpha * c[i+pntr[0]];
00146     for (j = pntr[i]; j < pntr[i+1]; j++)
00147       c[indx[j]] -= z * val[j];
00148       }
00149     else
00150       for (l = 0; l < n; l++)
00151     for (i = 0; i < m; i++) {
00152       z = alpha * _SpMatVal(c, ldc, i+pntr[0], l);
00153       for (j = pntr[i]; j < pntr[i+1]; j++)
00154         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00155     }
00156   }
00157 }
00158 
00159 
00160 static void
00161 CompCol_LowerUnitSolve_float(int m, int n, int unitd, const float *dv, 
00162     float alpha, 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   float z;
00167 
00168   // To make the compiler happy
00169   if (dv)
00170     ;
00171 
00172   if (unitd != 1) {
00173     cerr << "unitd != 1 not implemented" << endl;
00174     exit(1);
00175   }
00176 
00177   if (alpha == 0.0)
00178     return;
00179 
00180   CopyRectangularArray_float(m, n, b, ldb, &c[pntr[0]-1], ldc);
00181 
00182   c -= 1;
00183   val -= pntr[0];
00184   indx -= pntr[0];
00185 
00186   if (alpha == 1.0) {
00187     if (n == 1)
00188       for (i = 0; i < m; i++) {
00189     z = c[i+pntr[0]];
00190     for (j = pntr[i]; j < pntr[i+1]; j++)
00191       c[indx[j]] -= z * val[j];
00192       }
00193     else
00194       for (l = 0; l < n; l++)
00195     for (i = 0; i < m; i++) {
00196       z = _SpMatVal(c, ldc, i+pntr[0], l);
00197       for (j = pntr[i]; j < pntr[i+1]; j++)
00198         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00199     }
00200   } else {
00201     if (n == 1)
00202       for (i = 0; i < m; i++) {
00203     z = alpha * c[i+pntr[0]];
00204     for (j = pntr[i]; j < pntr[i+1]; j++)
00205       c[indx[j]] -= z * val[j];
00206       }
00207     else
00208       for (l = 0; l < n; l++)
00209     for (i = 0; i < m; i++) {
00210       z = alpha * _SpMatVal(c, ldc, i+pntr[0], l);
00211       for (j = pntr[i]; j < pntr[i+1]; j++)
00212         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00213     }
00214   }
00215 }
00216 
00217 
00218 static void
00219 CompCol_LowerDiagSolve_double(int m, int n, int unitd, const double *dv, double alpha,
00220                const double *val, const int *indx, const int *pntr,
00221                const double *b, int ldb, double *c, int ldc)
00222 {
00223   int i, j, l;
00224   double z;
00225 
00226   // To make the compiler happy
00227   if (dv)
00228     ;
00229 
00230   if (unitd != 1) {
00231     cerr << "unitd != 1 not implemented" << endl;
00232     exit(1);
00233   }
00234 
00235   if (alpha == 0.0)
00236     return;
00237 
00238   CopyRectangularArray_double(m, n, b, ldb, &c[pntr[0]-1], ldc);
00239 
00240   c -= 1;
00241   val -= pntr[0];
00242   indx -= pntr[0];
00243 
00244   if (alpha == 1.0) {
00245     if (n == 1)
00246       for (i = 0; i < m; i++) {
00247     z = c[i+pntr[0]] / val[pntr[i]];
00248     c[i+pntr[0]] = z;
00249     for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00250       c[indx[j]] -= z * val[j];
00251       }
00252     else
00253       for (l = 0; l < n; l++)
00254     for (i = 0; i < m; i++) {
00255       z = _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i]];
00256       _SpMatVal(c, ldc, i+pntr[0], l) = z;
00257       for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00258         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00259     }
00260   } else {
00261     if (n == 1)
00262       for (i = 0; i < m; i++) {
00263     z = alpha * c[i+pntr[0]] / val[pntr[i]];
00264     c[i+pntr[0]] = z;
00265     for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00266       c[indx[j]] -= z * val[j];
00267       }
00268     else
00269       for (l = 0; l < n; l++)
00270     for (i = 0; i < m; i++) {
00271       z = alpha * _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i]];
00272       _SpMatVal(c, ldc, i, l) = z;
00273       for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00274         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00275     }
00276   }
00277 }
00278 
00279 static void
00280 CompCol_LowerDiagSolve_float(int m, int n, int unitd, const float *dv, 
00281                 float alpha,
00282                const float *val, const int *indx, const int *pntr,
00283                const float *b, int ldb, float *c, int ldc)
00284 {
00285   int i, j, l;
00286   float z;
00287 
00288   // To make the compiler happy
00289   if (dv)
00290     ;
00291 
00292   if (unitd != 1) {
00293     cerr << "unitd != 1 not implemented" << endl;
00294     exit(1);
00295   }
00296 
00297   if (alpha == 0.0)
00298     return;
00299 
00300   CopyRectangularArray_float(m, n, b, ldb, &c[pntr[0]-1], ldc);
00301 
00302   c -= 1;
00303   val -= pntr[0];
00304   indx -= pntr[0];
00305 
00306   if (alpha == 1.0) {
00307     if (n == 1)
00308       for (i = 0; i < m; i++) {
00309     z = c[i+pntr[0]] / val[pntr[i]];
00310     c[i+pntr[0]] = z;
00311     for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00312       c[indx[j]] -= z * val[j];
00313       }
00314     else
00315       for (l = 0; l < n; l++)
00316     for (i = 0; i < m; i++) {
00317       z = _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i]];
00318       _SpMatVal(c, ldc, i+pntr[0], l) = z;
00319       for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00320         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00321     }
00322   } else {
00323     if (n == 1)
00324       for (i = 0; i < m; i++) {
00325     z = alpha * c[i+pntr[0]] / val[pntr[i]];
00326     c[i+pntr[0]] = z;
00327     for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00328       c[indx[j]] -= z * val[j];
00329       }
00330     else
00331       for (l = 0; l < n; l++)
00332     for (i = 0; i < m; i++) {
00333       z = alpha * _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i]];
00334       _SpMatVal(c, ldc, i, l) = z;
00335       for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00336         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00337     }
00338   }
00339 }
00340 
00341 
00342 static void
00343 CompCol_UpperUnitSolve_double(int m, int n, int unitd, const double *dv, 
00344                 double alpha,
00345                const double *val, const int *indx, const int *pntr,
00346                const double *b, int ldb, double *c, int ldc)
00347 {
00348   int i, j, l;
00349   double z;
00350 
00351   // To make the compiler happy
00352   if (dv)
00353     ;
00354 
00355   if (unitd != 1) {
00356     cerr << "unitd != 1 not implemented" << endl;
00357     exit(1);
00358   }
00359 
00360   if (alpha == 0.0)
00361     return;
00362 
00363   CopyRectangularArray_double(m, n, b, ldb, &c[pntr[0]-1], ldc);
00364 
00365   c -= 1;
00366   val -= pntr[0];  
00367   indx -= pntr[0];
00368 
00369   if (alpha == 1.0) {
00370     if (n == 1)
00371       for (i = m - 1; i >= 0; i--) {
00372     z = c[i+pntr[0]];
00373     for (j = pntr[i]; j < pntr[i+1]; j++)
00374       c[indx[j]] -= z * val[j];
00375       }
00376     else
00377       for (l = 0; l < n; l++)
00378     for (i = m - 1; i >= 0; i--) {
00379       z = _SpMatVal(c, ldc, i+pntr[0], l);
00380       for (j = pntr[i]; j < pntr[i+1]; j++)
00381         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00382     }
00383   } else {
00384     if (n == 1)
00385       for (i = m - 1; i >= 0; i--) {
00386     z = alpha * c[i+pntr[0]];
00387     for (j = pntr[i]; j < pntr[i+1]; j++)
00388       c[indx[j]] -= z * val[j];
00389       }
00390     else
00391       for (l = 0; l < n; l++)
00392     for (i = m - 1; i >= 0; i--) {
00393       z = alpha * _SpMatVal(c, ldc, i+pntr[0], l);
00394       for (j = pntr[i]; j < pntr[i+1]; j++)
00395         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00396     }
00397   }
00398 }
00399 
00400 
00401 static void
00402 CompCol_UpperUnitSolve_float(int m, int n, int unitd, const float *dv, 
00403                 float alpha,
00404                const float *val, const int *indx, const int *pntr,
00405                const float *b, int ldb, float *c, int ldc)
00406 {
00407   int i, j, l;
00408   float z;
00409 
00410   // To make the compiler happy
00411   if (dv)
00412     ;
00413 
00414   if (unitd != 1) {
00415     cerr << "unitd != 1 not implemented" << endl;
00416     exit(1);
00417   }
00418 
00419   if (alpha == 0.0)
00420     return;
00421 
00422   CopyRectangularArray_float(m, n, b, ldb, &c[pntr[0]-1], ldc);
00423 
00424   c -= 1;
00425   val -= pntr[0];  
00426   indx -= pntr[0];
00427 
00428   if (alpha == 1.0) {
00429     if (n == 1)
00430       for (i = m - 1; i >= 0; i--) {
00431     z = c[i+pntr[0]];
00432     for (j = pntr[i]; j < pntr[i+1]; j++)
00433       c[indx[j]] -= z * val[j];
00434       }
00435     else
00436       for (l = 0; l < n; l++)
00437     for (i = m - 1; i >= 0; i--) {
00438       z = _SpMatVal(c, ldc, i+pntr[0], l);
00439       for (j = pntr[i]; j < pntr[i+1]; j++)
00440         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00441     }
00442   } else {
00443     if (n == 1)
00444       for (i = m - 1; i >= 0; i--) {
00445     z = alpha * c[i+pntr[0]];
00446     for (j = pntr[i]; j < pntr[i+1]; j++)
00447       c[indx[j]] -= z * val[j];
00448       }
00449     else
00450       for (l = 0; l < n; l++)
00451     for (i = m - 1; i >= 0; i--) {
00452       z = alpha * _SpMatVal(c, ldc, i+pntr[0], l);
00453       for (j = pntr[i]; j < pntr[i+1]; j++)
00454         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00455     }
00456   }
00457 }
00458 
00459 
00460 static void
00461 CompCol_UpperDiagSolve_double(int m, int n, int unitd, const double *dv, 
00462             double alpha,
00463                const double *val, const int *indx, const int *pntr,
00464                const double *b, int ldb, double *c, int ldc)
00465 {
00466   int i, j, l;
00467   double z;
00468 
00469   // To make the compiler happy
00470   if (dv)
00471     ;
00472 
00473   if (unitd != 1) {
00474     cerr << "unitd != 1 not implemented" << endl;
00475     exit(1);
00476   }
00477 
00478   if (alpha == 0.0)
00479     return;
00480 
00481   CopyRectangularArray_double(m, n, b, ldb, &c[pntr[0]-1], ldc);
00482 
00483   c -= 1;
00484   val -= pntr[0];
00485   indx -= pntr[0];
00486 
00487   if (alpha == 1.0) {
00488     if (n == 1)
00489       for (i = m - 1; i >= 0; i--) {
00490     z = c[i+pntr[0]] / val[pntr[i+1]-1];
00491     c[i+pntr[0]] = z;
00492     for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00493       c[indx[j]] -= z * val[j];
00494       }
00495     else
00496       for (l = 0; l < n; l++)
00497     for (i = m - 1; i >= 0; i--) {
00498       z = _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i+1]-1];
00499       _SpMatVal(c, ldc, i+pntr[0], l) = z;
00500       for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00501         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00502     }
00503   } else {
00504     if (n == 1)
00505       for (i = m - 1; i >= 0; i--) {
00506     z = alpha * c[i+pntr[0]] / val[pntr[i+1]-1];
00507     c[i+pntr[0]] = z;
00508     for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00509       c[indx[j]] -= z * val[j];
00510       }
00511     else
00512       for (l = 0; l < n; l++)
00513     for (i = m - 1; i >= 0; i--) {
00514       z = alpha * _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i+1]-1];
00515       _SpMatVal(c, ldc, i+pntr[0], l) = z;
00516       for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00517         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00518     }
00519   }
00520 }
00521 
00522 
00523 static void
00524 CompCol_UpperDiagSolve_float(int m, int n, int unitd, const float *dv, 
00525                 float alpha,
00526                const float *val, const int *indx, const int *pntr,
00527                const float *b, int ldb, float *c, int ldc)
00528 {
00529   int i, j, l;
00530   float z;
00531 
00532   // To make the compiler happy
00533   if (dv)
00534     ;
00535 
00536   if (unitd != 1) {
00537     cerr << "unitd != 1 not implemented" << endl;
00538     exit(1);
00539   }
00540 
00541   if (alpha == 0.0)
00542     return;
00543 
00544   CopyRectangularArray_float(m, n, b, ldb, &c[pntr[0]-1], ldc);
00545 
00546   c -= 1;
00547   val -= pntr[0];
00548   indx -= pntr[0];
00549 
00550   if (alpha == 1.0) {
00551     if (n == 1)
00552       for (i = m - 1; i >= 0; i--) {
00553     z = c[i+pntr[0]] / val[pntr[i+1]-1];
00554     c[i+pntr[0]] = z;
00555     for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00556       c[indx[j]] -= z * val[j];
00557       }
00558     else
00559       for (l = 0; l < n; l++)
00560     for (i = m - 1; i >= 0; i--) {
00561       z = _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i+1]-1];
00562       _SpMatVal(c, ldc, i+pntr[0], l) = z;
00563       for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00564         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00565     }
00566   } else {
00567     if (n == 1)
00568       for (i = m - 1; i >= 0; i--) {
00569     z = alpha * c[i+pntr[0]] / val[pntr[i+1]-1];
00570     c[i+pntr[0]] = z;
00571     for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00572       c[indx[j]] -= z * val[j];
00573       }
00574     else
00575       for (l = 0; l < n; l++)
00576     for (i = m - 1; i >= 0; i--) {
00577       z = alpha * _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i+1]-1];
00578       _SpMatVal(c, ldc, i+pntr[0], l) = z;
00579       for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00580         _SpMatVal(c, ldc, indx[j], l) -= z * val[j];
00581     }
00582   }
00583 }
00584 
00585 
00586 static void
00587 CompRow_LowerUnitSolve_double(int m, int n, int unitd, const double *dv, 
00588                 double alpha,
00589                const double *val, const int *indx, const int *pntr,
00590                const double *b, int ldb, double *c, int ldc)
00591 {
00592   int i, j, l;
00593   double z;
00594 
00595   // To make the compiler happy
00596   if (dv)
00597     ;
00598 
00599   if (unitd != 1) {
00600     cerr << "unitd != 1 not implemented" << endl;
00601     exit(1);
00602   }
00603 
00604   if (alpha == 0.0)
00605     return;
00606 
00607   c -= 1;
00608   val -= pntr[0];
00609   indx -= pntr[0];
00610 
00611   if (alpha == 1.0) {
00612     if (n == 1)
00613       for (i = 0; i < m; i++) {
00614     z = 0;
00615     for (j = pntr[i]; j < pntr[i+1]; j++)
00616       z += c[indx[j]] * val[j];
00617     c[i+pntr[0]] = b[i] - z;
00618       }
00619     else
00620       for (l = 0; l < n; l++)
00621     for (i = 0; i < m; i++) {
00622       z = 0;
00623       for (j = pntr[i]; j < pntr[i+1]; j++)
00624         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00625       _SpMatVal(c, ldc, i+pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);
00626     }
00627   } else {
00628     if (n == 1)
00629       for (i = 0; i < m; i++) {
00630     z = 0;
00631     for (j = pntr[i]; j < pntr[i+1]; j++)
00632       z += c[indx[j]] * val[j];
00633     c[i+1] = alpha * (b[i] - z);
00634       }
00635     else
00636       for (l = 0; l < n; l++)
00637     for (i = 0; i < m; i++) {
00638       z = 0;
00639       for (j = pntr[i]; j < pntr[i+1]; j++)
00640         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00641       _SpMatVal(c, ldc, i+pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);
00642     }
00643   }
00644 }
00645 
00646 
00647 static void
00648 CompRow_LowerUnitSolve_float(int m, int n, int unitd, const float *dv, 
00649                 float alpha,
00650                const float *val, const int *indx, const int *pntr,
00651                const float *b, int ldb, float *c, int ldc)
00652 {
00653   int i, j, l;
00654   float z;
00655 
00656   // To make the compiler happy
00657   if (dv)
00658     ;
00659 
00660   if (unitd != 1) {
00661     cerr << "unitd != 1 not implemented" << endl;
00662     exit(1);
00663   }
00664 
00665   if (alpha == 0.0)
00666     return;
00667 
00668   c -= 1;
00669   val -= pntr[0];
00670   indx -= pntr[0];
00671 
00672   if (alpha == 1.0) {
00673     if (n == 1)
00674       for (i = 0; i < m; i++) {
00675     z = 0;
00676     for (j = pntr[i]; j < pntr[i+1]; j++)
00677       z += c[indx[j]] * val[j];
00678     c[i+pntr[0]] = b[i] - z;
00679       }
00680     else
00681       for (l = 0; l < n; l++)
00682     for (i = 0; i < m; i++) {
00683       z = 0;
00684       for (j = pntr[i]; j < pntr[i+1]; j++)
00685         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00686       _SpMatVal(c, ldc, i+pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);
00687     }
00688   } else {
00689     if (n == 1)
00690       for (i = 0; i < m; i++) {
00691     z = 0;
00692     for (j = pntr[i]; j < pntr[i+1]; j++)
00693       z += c[indx[j]] * val[j];
00694     c[i+1] = alpha * (b[i] - z);
00695       }
00696     else
00697       for (l = 0; l < n; l++)
00698     for (i = 0; i < m; i++) {
00699       z = 0;
00700       for (j = pntr[i]; j < pntr[i+1]; j++)
00701         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00702       _SpMatVal(c, ldc, i+pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);
00703     }
00704   }
00705 }
00706 
00707 
00708 static void
00709 CompRow_LowerDiagSolve_double(int m, int n, int unitd, const double *dv, 
00710                 double alpha,
00711                const double *val, const int *indx, const int *pntr,
00712                const double *b, int ldb, double *c, int ldc)
00713 {
00714   int i, j, l;
00715   double z;
00716 
00717   // To make the compiler happy
00718   if (dv)
00719     ;
00720 
00721   if (unitd != 1) {
00722     cerr << "unitd != 1 not implemented" << endl;
00723     exit(1);
00724   }
00725 
00726   if (alpha == 0.0)
00727     return;
00728 
00729   c -= 1;
00730   val -= pntr[0];
00731   indx -= pntr[0];
00732 
00733   if (alpha == 1.0) {
00734     if (n == 1)
00735       for (i = 0; i < m; i++) {
00736     z = 0;
00737     for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00738       z += c[indx[j]] * val[j];
00739     c[i+pntr[0]] = (b[i] - z) / val[pntr[i+1]-1];
00740       }
00741     else
00742       for (l = 0; l < n; l++)
00743     for (i = 0; i < m; i++) {
00744       z = 0;
00745       for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00746         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00747       _SpMatVal(c, ldc, i+pntr[0], l) = 
00748         (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i+1]-1];
00749     }
00750   } else {
00751     if (n == 1)
00752       for (i = 0; i < m; i++) {
00753     z = 0;
00754     for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00755       z += c[indx[j]] * val[j];
00756     c[i+pntr[0]] = alpha * (b[i] - z) / val[pntr[i+1]-1];
00757       }
00758     else
00759       for (l = 0; l < n; l++)
00760     for (i = 0; i < m; i++) {
00761       z = 0;
00762       for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00763         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00764       _SpMatVal(c, ldc, i+pntr[0], l) = 
00765         alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i+1]-1];
00766     }
00767   }
00768 }
00769 
00770 
00771 static void
00772 CompRow_LowerDiagSolve_float(int m, int n, int unitd, const float *dv, 
00773                 float alpha,
00774                const float *val, const int *indx, const int *pntr,
00775                const float *b, int ldb, float *c, int ldc)
00776 {
00777   int i, j, l;
00778   float z;
00779 
00780   // To make the compiler happy
00781   if (dv)
00782     ;
00783 
00784   if (unitd != 1) {
00785     cerr << "unitd != 1 not implemented" << endl;
00786     exit(1);
00787   }
00788 
00789   if (alpha == 0.0)
00790     return;
00791 
00792   c -= 1;
00793   val -= pntr[0];
00794   indx -= pntr[0];
00795 
00796   if (alpha == 1.0) {
00797     if (n == 1)
00798       for (i = 0; i < m; i++) {
00799     z = 0;
00800     for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00801       z += c[indx[j]] * val[j];
00802     c[i+pntr[0]] = (b[i] - z) / val[pntr[i+1]-1];
00803       }
00804     else
00805       for (l = 0; l < n; l++)
00806     for (i = 0; i < m; i++) {
00807       z = 0;
00808       for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00809         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00810       _SpMatVal(c, ldc, i+pntr[0], l) = 
00811         (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i+1]-1];
00812     }
00813   } else {
00814     if (n == 1)
00815       for (i = 0; i < m; i++) {
00816     z = 0;
00817     for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00818       z += c[indx[j]] * val[j];
00819     c[i+pntr[0]] = alpha * (b[i] - z) / val[pntr[i+1]-1];
00820       }
00821     else
00822       for (l = 0; l < n; l++)
00823     for (i = 0; i < m; i++) {
00824       z = 0;
00825       for (j = pntr[i]; j < pntr[i+1] - 1; j++)
00826         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00827       _SpMatVal(c, ldc, i+pntr[0], l) = 
00828         alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i+1]-1];
00829     }
00830   }
00831 }
00832 
00833 
00834 static void
00835 CompRow_UpperUnitSolve_double(int m, int n, int unitd, const double *dv, double alpha,
00836                const double *val, const int *indx, const int *pntr,
00837                const double *b, int ldb, double *c, int ldc)
00838 
00839 {
00840   int i, j, l;
00841   double z;
00842 
00843   // To make the compiler happy
00844   if (dv)
00845     ;
00846 
00847   if (unitd != 1) {
00848     cerr << "unitd != 1 not implemented" << endl;
00849     exit(1);
00850   }
00851 
00852   if (alpha == 0.0)
00853     return;
00854 
00855   c -= 1;
00856   val -= pntr[0];
00857   indx -= pntr[0];
00858 
00859   if (alpha == 1.0) {
00860     if (n == 1)
00861       for (i = m - 1; i >= 0; i--) {
00862     z = 0;
00863     for (j = pntr[i]; j < pntr[i+1]; j++)
00864       z += c[indx[j]] * val[j];
00865     c[i+pntr[0]] = b[i] - z;
00866       }
00867     else
00868       for (l = 0; l < n; l++)
00869     for (i = m - 1; i >= 0; i--) {
00870       z = 0;
00871       for (j = pntr[i]; j < pntr[i+1]; j++)
00872         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00873       _SpMatVal(c, ldc, i+pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);
00874     }
00875   } else {
00876     if (n == 1)
00877       for (i = m - 1; i >= 0; i--) {
00878     z = 0;
00879     for (j = pntr[i]; j < pntr[i+1]; j++)
00880       z += c[indx[j]] * val[j];
00881     c[i+pntr[0]] = alpha * (b[i] - z);
00882       }
00883     else
00884       for (l = 0; l < n; l++)
00885     for (i = m - 1; i >= 0; i--) {
00886       z = 0;
00887       for (j = pntr[i]; j < pntr[i+1]; j++)
00888         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00889       _SpMatVal(c, ldc, i+pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);
00890     }
00891   }
00892 }
00893 
00894 
00895 static void
00896 CompRow_UpperUnitSolve_float(int m, int n, int unitd, const float *dv, float alpha,
00897                const float *val, const int *indx, const int *pntr,
00898                const float *b, int ldb, float *c, int ldc)
00899 
00900 {
00901   int i, j, l;
00902   float z;
00903 
00904   // To make the compiler happy
00905   if (dv)
00906     ;
00907 
00908   if (unitd != 1) {
00909     cerr << "unitd != 1 not implemented" << endl;
00910     exit(1);
00911   }
00912 
00913   if (alpha == 0.0)
00914     return;
00915 
00916   c -= 1;
00917   val -= pntr[0];
00918   indx -= pntr[0];
00919 
00920   if (alpha == 1.0) {
00921     if (n == 1)
00922       for (i = m - 1; i >= 0; i--) {
00923     z = 0;
00924     for (j = pntr[i]; j < pntr[i+1]; j++)
00925       z += c[indx[j]] * val[j];
00926     c[i+pntr[0]] = b[i] - z;
00927       }
00928     else
00929       for (l = 0; l < n; l++)
00930     for (i = m - 1; i >= 0; i--) {
00931       z = 0;
00932       for (j = pntr[i]; j < pntr[i+1]; j++)
00933         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00934       _SpMatVal(c, ldc, i+pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);
00935     }
00936   } else {
00937     if (n == 1)
00938       for (i = m - 1; i >= 0; i--) {
00939     z = 0;
00940     for (j = pntr[i]; j < pntr[i+1]; j++)
00941       z += c[indx[j]] * val[j];
00942     c[i+pntr[0]] = alpha * (b[i] - z);
00943       }
00944     else
00945       for (l = 0; l < n; l++)
00946     for (i = m - 1; i >= 0; i--) {
00947       z = 0;
00948       for (j = pntr[i]; j < pntr[i+1]; j++)
00949         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00950       _SpMatVal(c, ldc, i+pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);
00951     }
00952   }
00953 }
00954 
00955 
00956 static void
00957 CompRow_UpperDiagSolve_double(int m, int n, int unitd, const double *dv, double alpha,
00958                const double *val, const int *indx, const int *pntr,
00959                const double *b, int ldb, double *c, int ldc)
00960 {
00961   int i, j, l;
00962   double z;
00963 
00964   // To make the compiler happy
00965   if (dv)
00966     ;
00967 
00968   if (unitd != 1) {
00969     cerr << "unitd != 1 not implemented" << endl;
00970     exit(1);
00971   }
00972 
00973   if (alpha == 0.0)
00974     return;
00975 
00976   c -= 1;
00977   val -= pntr[0];
00978   indx -= pntr[0];
00979   
00980   if (alpha == 1.0) {
00981     if (n == 1)
00982       for (i = m - 1; i >= 0; i--) {    
00983     z = 0;
00984     for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00985       z += c[indx[j]] * val[j];
00986     c[i+pntr[0]] = (b[i] - z) / val[pntr[i]];
00987       }
00988     else
00989       for (l = 0; l < n; l++)
00990     for (i = m - 1; i >= 0; i--) {  
00991       z = 0;
00992       for (j = pntr[i] + 1; j < pntr[i+1]; j++)
00993         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
00994       _SpMatVal(c, ldc, i+pntr[0], l) = 
00995         (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i]];
00996     }
00997   } else {
00998     if (n == 1)
00999       for (i = m - 1; i >= 0; i--) {    
01000     z = 0;
01001     for (j = pntr[i] + 1; j < pntr[i+1]; j++)
01002       z += c[indx[j]] * val[j];
01003     c[i+pntr[0]] = alpha * (b[i] - z) / val[pntr[i]];
01004       }
01005     else
01006       for (l = 0; l < n; l++)
01007     for (i = m - 1; i >= 0; i--) {  
01008       z = 0;
01009       for (j = pntr[i] + 1; j < pntr[i+1]; j++)
01010         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
01011       _SpMatVal(c, ldc, i+pntr[0], l) = 
01012         alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i]];
01013     }
01014   }
01015 }
01016 
01017 
01018 static void
01019 CompRow_UpperDiagSolve_float(int m, int n, int unitd, const float *dv, float alpha,
01020                const float *val, const int *indx, const int *pntr,
01021                const float *b, int ldb, float *c, int ldc)
01022 {
01023   int i, j, l;
01024   float z;
01025 
01026   // To make the compiler happy
01027   if (dv)
01028     ;
01029 
01030   if (unitd != 1) {
01031     cerr << "unitd != 1 not implemented" << endl;
01032     exit(1);
01033   }
01034 
01035   if (alpha == 0.0)
01036     return;
01037 
01038   c -= 1;
01039   val -= pntr[0];
01040   indx -= pntr[0];
01041   
01042   if (alpha == 1.0) {
01043     if (n == 1)
01044       for (i = m - 1; i >= 0; i--) {    
01045     z = 0;
01046     for (j = pntr[i] + 1; j < pntr[i+1]; j++)
01047       z += c[indx[j]] * val[j];
01048     c[i+pntr[0]] = (b[i] - z) / val[pntr[i]];
01049       }
01050     else
01051       for (l = 0; l < n; l++)
01052     for (i = m - 1; i >= 0; i--) {  
01053       z = 0;
01054       for (j = pntr[i] + 1; j < pntr[i+1]; j++)
01055         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
01056       _SpMatVal(c, ldc, i+pntr[0], l) = 
01057         (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i]];
01058     }
01059   } else {
01060     if (n == 1)
01061       for (i = m - 1; i >= 0; i--) {    
01062     z = 0;
01063     for (j = pntr[i] + 1; j < pntr[i+1]; j++)
01064       z += c[indx[j]] * val[j];
01065     c[i+pntr[0]] = alpha * (b[i] - z) / val[pntr[i]];
01066       }
01067     else
01068       for (l = 0; l < n; l++)
01069     for (i = m - 1; i >= 0; i--) {  
01070       z = 0;
01071       for (j = pntr[i] + 1; j < pntr[i+1]; j++)
01072         z += _SpMatVal(c, ldc, indx[j], l) * val[j];
01073       _SpMatVal(c, ldc, i+pntr[0], l) = 
01074         alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i]];
01075     }
01076   }
01077 }
01078 
01079 
01080 
01081 
01082 
01083 
01084 
01085 /*
01086  * C <- alpha D A^{-1} B + beta C
01087  * C <- alpha A^{-1} D B + beta C
01088  */
01089 IFPACK_DEPRECATED void F77NAME(scscsm)
01090   (const int &transa, const int &m, const int &n,
01091    const int &unitd, const float *dv, const float &alpha,
01092    const int descra[], const float *val,
01093    const int *indx, const int *pntr, const float *b, int &ldb,
01094    const float &beta, float *c, const int &ldc,
01095    float *work, const int &lwork)
01096 {
01097   if (descra[0] != 0) {
01098     cerr << "Must have general matrix" << endl;
01099     exit(1);
01100   }
01101 
01102   if (beta != 0.0) {
01103     cerr << "beta != 0 not implemented" << endl;
01104     exit(1);
01105   }
01106 
01107   // To make the compiler happy
01108   if (work && lwork)
01109     ;
01110 
01111   if (transa == 0) {
01112     if (descra[1] == 1) {
01113       if (descra[2] == 0)
01114     CompCol_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01115                    b, ldb, c, ldc);
01116       else if (descra[2] == 1)
01117     CompCol_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01118                    b, ldb, c, ldc);
01119       else {
01120     cerr << "Bad descra[2]" << endl;
01121     exit(1);
01122       }
01123     } else if (descra[1] == 2) {
01124       if (descra[2] == 0)
01125     CompCol_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01126                    b, ldb, c, ldc);
01127       else if (descra[2] == 1)
01128     CompCol_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01129                    b, ldb, c, ldc);
01130       else {
01131     cerr << "Bad descra[2]" << endl;
01132     exit(1);
01133       }
01134     }
01135   } else if (transa == 1 || transa == 2) {
01136     if (descra[1] == 1) {
01137       if (descra[2] == 0)
01138     CompRow_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01139                    b, ldb, c, ldc);
01140       else if (descra[2] == 1)
01141     CompRow_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01142                    b, ldb, c, ldc);
01143       else {
01144     cerr << "Bad descra[2]" << endl;
01145     exit(1);
01146       }
01147     } else if (descra[1] == 2) {
01148       if (descra[2] == 0)
01149     CompRow_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01150                    b, ldb, c, ldc);
01151       else if (descra[2] == 1)
01152     CompRow_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01153                    b, ldb, c, ldc);
01154       else {
01155     cerr << "Bad descra[2]" << endl;
01156     exit(1);
01157       }
01158     } else {
01159       cerr << "Bad descra[1]" << endl;
01160       exit(1);
01161     }
01162   } else {
01163     cerr << "Bad transa" << endl;
01164     exit(1);
01165   }
01166 }
01167 
01168 
01169 IFPACK_DEPRECATED void F77NAME(scsrsm)
01170   (const int &transa, const int &m, const int &n,
01171    const int &unitd, const float *dv, const float &alpha,
01172    const int descra[], const float *val,
01173    const int *indx, const int *pntr, const float *b, int &ldb,
01174    const float &beta, float *c, const int &ldc,
01175    float *work, const int &lwork)
01176 {
01177   if (descra[0] != 0) {
01178     cerr << "Must have general matrix" << endl;
01179     exit(1);
01180   }
01181 
01182   if (beta != 0.0) {
01183     cerr << "beta != 0 not implemented" << endl;
01184     exit(1);
01185   }
01186 
01187   // To make the compiler happy
01188   if (work && lwork)
01189     ;
01190 
01191   if (transa == 0) {
01192     if (descra[1] == 1) {
01193       if (descra[2] == 0)
01194     CompRow_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01195                    b, ldb, c, ldc);
01196       else if (descra[2] == 1)
01197     CompRow_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01198                    b, ldb, c, ldc);
01199       else {
01200     cerr << "Bad descra[2]" << endl;
01201     exit(1);
01202       }
01203     } else if (descra[1] == 2) {
01204       if (descra[2] == 0)
01205     CompRow_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01206                    b, ldb, c, ldc);
01207       else if (descra[2] == 1)
01208     CompRow_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01209                    b, ldb, c, ldc);
01210       else {
01211     cerr << "Bad descra[2]" << endl;
01212     exit(1);
01213       }
01214     }
01215   } else if (transa == 1 || transa == 2) {
01216     if (descra[1] == 1) {
01217       if (descra[2] == 0)
01218     CompCol_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01219                    b, ldb, c, ldc);
01220       else if (descra[2] == 1)
01221     CompCol_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01222                    b, ldb, c, ldc);
01223       else {
01224     cerr << "Bad descra[2]" << endl;
01225     exit(1);
01226       }
01227     } else if (descra[1] == 2) {
01228       if (descra[2] == 0)
01229     CompCol_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01230                    b, ldb, c, ldc);
01231       else if (descra[2] == 1)
01232     CompCol_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
01233                    b, ldb, c, ldc);
01234       else {
01235     cerr << "Bad descra[2]" << endl;
01236     exit(1);
01237       }
01238     } else {
01239       cerr << "Bad descra[1]" << endl;
01240       exit(1);
01241     }
01242   } else {
01243     cerr << "Bad transa" << endl;
01244     exit(1);
01245   }
01246 }
01247 
01248 
01249 IFPACK_DEPRECATED void F77NAME(dcscsm)
01250   (const int &transa, const int &m, const int &n,
01251    const int &unitd, const double *dv, const double &alpha,
01252    const int descra[], const double *val,
01253    const int *indx, const int *pntr, const double *b, int &ldb,
01254    const double &beta, double *c, const int &ldc,
01255    double *work, const int &lwork)
01256 {
01257   if (descra[0] != 0) {
01258     cerr << "Must have general matrix" << endl;
01259     exit(1);
01260   }
01261 
01262   if (beta != 0.0) {
01263     cerr << "beta != 0 not implemented" << endl;
01264     exit(1);
01265   }
01266 
01267   // To make the compiler happy
01268   if (work && lwork)
01269     ;
01270 
01271   if (transa == 0) {
01272     if (descra[1] == 1) {
01273       if (descra[2] == 0)
01274     CompCol_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01275                    b, ldb, c, ldc);
01276       else if (descra[2] == 1)
01277     CompCol_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01278                    b, ldb, c, ldc);
01279       else {
01280     cerr << "Bad descra[2]" << endl;
01281     exit(1);
01282       }
01283     } else if (descra[1] == 2) {
01284       if (descra[2] == 0)
01285     CompCol_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01286                    b, ldb, c, ldc);
01287       else if (descra[2] == 1)
01288     CompCol_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01289                    b, ldb, c, ldc);
01290       else {
01291     cerr << "Bad descra[2]" << endl;
01292     exit(1);
01293       }
01294     }
01295   } else if (transa == 1 || transa == 2) {
01296     if (descra[1] == 1) {
01297       if (descra[2] == 0)
01298     CompRow_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01299                    b, ldb, c, ldc);
01300       else if (descra[2] == 1)
01301     CompRow_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01302                    b, ldb, c, ldc);
01303       else {
01304     cerr << "Bad descra[2]" << endl;
01305     exit(1);
01306       }
01307     } else if (descra[1] == 2) {
01308       if (descra[2] == 0)
01309     CompRow_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01310                    b, ldb, c, ldc);
01311       else if (descra[2] == 1)
01312     CompRow_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01313                    b, ldb, c, ldc);
01314       else {
01315     cerr << "Bad descra[2]" << endl;
01316     exit(1);
01317       }
01318     } else {
01319       cerr << "Bad descra[1]" << endl;
01320       exit(1);
01321     }
01322   } else {
01323     cerr << "Bad transa" << endl;
01324     exit(1);
01325   }
01326 }
01327 
01328 
01329 
01330 
01331 IFPACK_DEPRECATED void F77NAME(dcsrsm)
01332   (const int &transa, const int &m, const int &n,
01333    const int &unitd, const double *dv, const double &alpha,
01334    const int descra[], const double *val,
01335    const int *indx, const int *pntr, const double *b, int &ldb,
01336    const double &beta, double *c, const int &ldc,
01337    double *work, const int &lwork)
01338 {
01339   if (descra[0] != 0) {
01340     cerr << "Must have general matrix" << endl;
01341     exit(1);
01342   }
01343 
01344   if (beta != 0.0) {
01345     cerr << "beta != 0 not implemented" << endl;
01346     exit(1);
01347   }
01348 
01349   // To make the compiler happy
01350   if (work && lwork)
01351     ;
01352 
01353   if (transa == 0) {
01354     if (descra[1] == 1) {
01355       if (descra[2] == 0)
01356     CompRow_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01357                    b, ldb, c, ldc);
01358       else if (descra[2] == 1)
01359     CompRow_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01360                    b, ldb, c, ldc);
01361       else {
01362     cerr << "Bad descra[2]" << endl;
01363     exit(1);
01364       }
01365     } else if (descra[1] == 2) {
01366       if (descra[2] == 0)
01367     CompRow_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01368                    b, ldb, c, ldc);
01369       else if (descra[2] == 1)
01370     CompRow_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01371                    b, ldb, c, ldc);
01372       else {
01373     cerr << "Bad descra[2]" << endl;
01374     exit(1);
01375       }
01376     }
01377   } else if (transa == 1 || transa == 2) {
01378     if (descra[1] == 1) {
01379       if (descra[2] == 0)
01380     CompCol_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01381                    b, ldb, c, ldc);
01382       else if (descra[2] == 1)
01383     CompCol_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01384                    b, ldb, c, ldc);
01385       else {
01386     cerr << "Bad descra[2]" << endl;
01387     exit(1);
01388       }
01389     } else if (descra[1] == 2) {
01390       if (descra[2] == 0)
01391     CompCol_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01392                    b, ldb, c, ldc);
01393       else if (descra[2] == 1)
01394     CompCol_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
01395                    b, ldb, c, ldc);
01396       else {
01397     cerr << "Bad descra[2]" << endl;
01398     exit(1);
01399       }
01400     } else {
01401       cerr << "Bad descra[1]" << endl;
01402       exit(1);
01403     }
01404   } else {
01405     cerr << "Bad transa" << endl;
01406     exit(1);
01407   }
01408 }
01409 
01410 
01411 
01412 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines