dcsrsm_c.c

用于求解大型稀疏线性方程组Ax=b的数值计算库.

/*-------------------------------------------------------|
|  NIST SPARSE BLAS v. 0.9 (Sat Jul 6 14:27:21 EDT 1996) |
|                                                        |
|  Authors:                                              |
|     Karin A. Remington and Roldan Pozo                 |
|     National Institute of Standards and Technology     |
|                                                        |
|  Based on the interface standard proposed in:          | 
|   "A Revised Proposal for a Sparse BLAS Toolkit" by    |
|    S. Carney and K. Wu -- University of Minnesota      |
|    M. Heroux and G. Li -- Cray Research                |  
|    R. Pozo and K.A. Remington -- NIST                  |
|                                                        |
|  Contact:                                              |
|     Karin A. Remington, email: [e-mail]     |
--------------------------------------------------------*/

#include <stdlib.h>
#include <stdio.h>
#include "spblas.h"
#include "dcscmtsl.h"
#include "dcscvtsl.h"
#include "dcsrmtsl.h"
#include "dcsrvtsl.h"


/* Sparse BLAS Toolkit interface routine: */
void  dcsrsm(
             const int transa, const int m, const int n, 
             const int unitd, const double dv[], 
             const double alpha, const int descra[], const double val[],
             const int indx[], const int pntrb[], const int pntre[],
             const double b[], const int ldb,
             const double beta, double c[], const int ldc,
             double work[], const int lwork)
{
/* ------------ begin interface description ------------
   Toolkit interface:
   dcsrsm -- compressed sparse row format triangular solve
  
   C <- alpha D inv(A) B + beta C    C <- alpha D inv(A') B + beta C
   C <- alpha inv(A) D B + beta C    C <- alpha inv(A') D B + beta C
   
                                      ( ' indicates matrix transpose)
  
   Arguments:
  
   int transa	Indicates how to operate with the sparse matrix
  		0 : operate with matrix
  		1 : operate with transpose matrix
  
   int m	Number of rows in matrix c
  
   int n	Number of columns in matrix c
  
   int unitd	Type of scaling:
                        1 : Identity matrix (argument dv[] is ignored)
                        2 : Scale on left (row scaling)
                        3 : Scale on right (column scaling)
  
   double alpha	Scalar parameter
  
   double beta 	Scalar parameter
  
   int descra[]	Descriptor argument.  Nine element integer array
  		descra[0] matrix structure
  			0 : general
  			1 : symmetric
  			2 : Hermitian
  			3 : Triangular
  			4 : Skew(Anti-Symmetric
  			5 : Diagonal
  		descra[1] upper/lower triangular indicator
  			1 : lower
  			2 : upper
  		descra[2] main diagonal type
  			0 : non-unit
  			1 : unit
  		descra[3] Array base 
  			0 : C/C++ compatible
  			1 : Fortran compatible
  		descra[4] repeated indices?
  			0 : unknown
  			1 : no repeated indices
  
  
   double *val	scalar array of length nnz containing matrix entries
  
   int *indx    integer array of length nnz containing column indices
  
   int *pntrb   integer array of length k such that pntrb(j)-pntrb(1)
                points to location in val of the first nonzero element in row j
  
   int *pntre   integer array of length k such that pntre(j)-pntrb(1)
                points to location in val of the last nonzero element in row j
  
   double *b	rectangular array with first dimension ldb
  
   double *c	rectangular array with first dimension ldc
  
   double *work	scratch array of length lwork.  
                lwork should be at least (n*m)
  
   ------------ end interface description --------------*/
int ind_base = descra[3];

  if (lwork < m*n ){
    printf("Insufficient work space for dcsrsm.\n");
    printf("   lwork must be at least (n*m) = %d \n",m*n);
    return;
  }

if (alpha == 0.0) {
   ScaleArray_double(m, n, c, ldc, beta);
   return;
} 
switch (descra[0]) {
case 3:  /* Matrix MUST be triangular, of course */
  switch (transa) {
  case 0:
    switch  ( 10*descra[1]+descra[2] ) {
    case 10:  /* Lower triangular, non-unit diagonal */
       switch (n) {
       case 1: /* Vec Mult */
         if (alpha == 1) {
           if (beta == 1) {
             switch (unitd) {
             case 1: /* No scaling */
                CSR_VecTriangSlvLD_CABC_double(m,  val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_VecTriangSlvLD_CDABC_double(m,  dv, val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_VecTriangSlvLD_CADBC_double(m,  dv, val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             default:
                printf("Invalid argument unitd in dcsrsm. Use 1-3. \n");
                break;
             } /* end switch on unitd */
#if (0)
           } else if (beta == -1) {
             switch (unitd) {
             case 1: /* No scaling */
                CSR_VecTriangSlvLD_CABmC_double(m,  val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_VecTriangSlvLD_CDABmC_double(m,  dv, val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_VecTriangSlvLD_CADBmC_double(m,  dv, val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             default:
                printf("Invalid argument unitd in dcsrsm. Use 1-3. \n");
                break;
             } /* end switch on unitd */
#endif
           } else if (beta == 0) {
             switch (unitd) {
             case 1: /* No scaling */
                CSR_VecTriangSlvLD_CAB_double(m,  val, indx, pntrb, pntre, b, 
                                     c, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_VecTriangSlvLD_CDAB_double(m,  dv, val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_VecTriangSlvLD_CADB_double(m,  dv, val, indx, pntrb, pntre, b, 
                                     c, ind_base);
                break;
             default:
                printf("Invalid argument unitd in dcsrsm. Use 1-3. \n");
                break;
             } /* end switch on unitd */
           } else  { /*  beta is general nonzero */
             switch (unitd) {
             case 1: /* No scaling */
                CSR_VecTriangSlvLD_CABbC_double(m,  val, indx, pntrb, pntre, b, 
                                     beta, c, work, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_VecTriangSlvLD_CDABbC_double(m,  dv, val, indx, pntrb, pntre, b, 
                                     beta, c, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_VecTriangSlvLD_CADBbC_double(m,  dv, val, indx, pntrb, pntre, b, 
                                     beta, c, work, ind_base);
                break;
             default:
                printf("Invalid argument unitd in dcsrsm. Use 1-3. \n");
                break;
             } /* end switch on unitd */
           }
         } else { /* alpha is general nonzero */
           if (beta == 1) {
             switch (unitd) {
             case 1: /* No scaling */
                CSR_VecTriangSlvLD_CaABC_double(m,  alpha, val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_VecTriangSlvLD_CaDABC_double(m,  dv, alpha, val, indx, pntrb, 
                                     pntre, b, c, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_VecTriangSlvLD_CaADBC_double(m,  dv, alpha, val, indx, pntrb, 
                                     pntre, b, c, work, ind_base);
                break;
             default:
                printf("Invalid argument unitd in dcsrsm. Use 1-3. \n");
                break;
             } /* end switch on unitd */
#if (0)
           } else if (beta == -1) {
             switch (unitd) {
             case 1: /* No scaling */
                CSR_VecTriangSlvLD_CaABmC_double(m,  alpha, val, indx, pntrb, pntre, b, 
                                     c, work, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_VecTriangSlvLD_CaDABmC_double(m,  dv, alpha, val, indx, pntrb, 
                                     pntre, b, c, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_VecTriangSlvLD_CaADBmC_double(m,  dv, alpha, val, indx, pntrb, 
                                     pntre, b, c, work, ind_base);
                break;
             default:
                printf("Invalid argument unitd in dcsrsm. Use 1-3. \n");
                break;
             } /* end switch on unitd */
#endif
           } else if (beta == 0) {
             switch (unitd) {
             case 1: /* No scaling */
                CSR_VecTriangSlvLD_CaAB_double(m,  alpha, val, indx, pntrb, pntre, b, 
                                     c, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_VecTriangSlvLD_CaDAB_double(m,  dv, alpha, val, indx, pntrb, 
                                     pntre, b, c, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_VecTriangSlvLD_CaADB_double(m,  dv, alpha, val, indx, pntrb, 
                                     pntre, b, c, ind_base);
                break;
             default:
                printf("Invalid argument unitd in dcsrsm. Use 1-3. \n");
                break;
             } /* end switch on unitd */
           } else { /*  beta is general nonzero */
             switch (unitd) {
             case 1: /* No scaling */
                CSR_VecTriangSlvLD_CaABbC_double(m,  alpha, val, indx, pntrb, pntre, 
                                     b, beta, c, work, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_VecTriangSlvLD_CaDABbC_double(m,  dv, alpha, val, indx, pntrb, 
                                     pntre, b, beta, c, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_VecTriangSlvLD_CaADBbC_double(m,  dv, alpha, val, indx, pntrb, 
                                     pntre, b, beta, c, work, ind_base);
                break;
             default:
                printf("Invalid argument unitd in dcsrsm. Use 1-3. \n");
                break;
             } /* end switch on unitd */
           }
         }
         break;
       default: /* Mat Mult */
         if (alpha == 1) {
           if (beta == 1) {
             switch (unitd) {
             case 1: /* No scaling */
                CSR_MatTriangSlvLD_CABC_double(m, n,  val, indx, pntrb, pntre, b, 
                                     ldb, c, ldc, work, ind_base);
                break;
             case 2: /* Left scaling */
                CSR_MatTriangSlvLD_CDABC_double(m, n,  dv, val, indx, pntrb, pntre, 
                                     b, ldb, c, ldc, work, ind_base);
                break;
             case 3: /* Right scaling */
                CSR_MatTriangSlvLD_CADBC_double(m, n,  dv, val, indx, pntrb, pntre, 
                                     b, ldb, c, ldc, work, ind_base);