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