|
IFPACK Development
|
|
IFPACK Development
|
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
1.7.4
