dvbrmtsl.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]     |
--------------------------------------------------------*/


/* Created:  Sat Jul 6 14:38:13 EDT 1996 */

#include "dvbrmtsl.h"



void VBR_MatTriangSlvLU_CAB_double(
                 const int mb, const int n, 
                  
                 
                 const double *val, const int *indx,
                 const int *bindx, const int *rpntr,
                 const int *cpntr, const int *bpntrb, const int *bpntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc,
                  const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,mm,bs,br,cs,cr;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=rpntr[mb]-rpntr[0];
  double t;

  val-=ind_base;
  b-=ind_base;
  c-=ind_base;
  indx-=ind_base;
  bindx-=ind_base;
  cpntr-=ind_base;



  
  
  mm=0; 
  ds=0;
  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    ds += mm*mm;                
    mm = rpntr[i+1]-rpntr[i];   
    pb = &b[bs];                
    pc = &c[bs];                

    

    for (l=0;l!=n;l++){                                      
      for (j=0;j!=mm;j++) pc[j] = pb[j];
      pc+=ldc;pb+=ldb;                                       
    }                                                        
  }

  

  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    mm = rpntr[i+1]-rpntr[i];   

    for (j=jb;j!=je;j++) {      
      index = bindx[j];             
      if ( index == i + ind_base ){ 
        continue;                   
      } else {
        pval = &val[indx[j]];       
        ptmp = pval;                                         
                                 
        cs = cpntr[index];          
        cr = cpntr[index+1] - cs;   
        pb = &b[bs];                
        pc = &c[bs];                
        cj = &c[cs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=cr;ii++) {      
            t = cj[ii];
            for (jj=0;jj!=mm;jj++) {        
              pc[jj] -= (*pval++) * t;       
            }
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }
  }

  
  

        



}


void VBR_MatTriangSlvUU_CAB_double(
                 const int mb, const int n,
                  
                 
                 const double *val, const int *indx,
                 const int *bindx, const int *rpntr,
                 const int *cpntr, const int *bpntrb, const int *bpntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc,
                  const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,mm,bs,br,cs,cr;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=rpntr[mb]-rpntr[0];
  double t;

  val-=ind_base;
  b-=ind_base;
  c-=ind_base;
  indx-=ind_base;
  bindx-=ind_base;
  cpntr-=ind_base;



  
  
  mm=0; 
  ds=0;
  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    ds += mm*mm;                
    mm = rpntr[i+1]-rpntr[i];   
    pb = &b[bs];                
    pc = &c[bs];                

    

    for (l=0;l!=n;l++){                                      
      for (j=0;j!=mm;j++) pc[j] = pb[j];
      pc+=ldc;pb+=ldb;                                       
    }                                                        
  }

  

  for (i=mb-1;i!=-1;i--) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    mm = rpntr[i+1]-rpntr[i];   
    pb = &b[bs];                
    pc = &c[bs];                

    for (j=jb;j!=je;j++) {      
      index = bindx[j];             
      if ( index == i + ind_base ){ 
        continue;                   
      } else {
        pval = &val[indx[j]];       
        ptmp = pval;                                         
                                 
        cs = cpntr[index];          
        cr = cpntr[index+1] - cs;   
        pb = &b[bs];                
        pc = &c[bs];                
        cj = &c[cs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=cr;ii++) {      
            t = cj[ii];
            for (jj=0;jj!=mm;jj++) {        
              pc[jj] -= (*pval++) * t;       
            }
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }
  }

  
  

        




}

void VBR_MatTriangSlvLU_CATB_double(
                 const int mb, const int n,
                  
                 
                 const double *val, const int *indx,
                 const int *bindx, const int *rpntr,
                 const int *cpntr, const int *bpntrb, const int *bpntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc,
                  const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,mm,bs,br,cs,cr;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=rpntr[mb]-rpntr[0];
  double t;

  val-=ind_base;
  b-=ind_base;
  c-=ind_base;
  indx-=ind_base;
  bindx-=ind_base;
  cpntr-=ind_base;



  
  
  mm=0; 
  ds=0;
  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    ds += mm*mm;                
    mm = rpntr[i+1]-rpntr[i];   
    pb = &b[bs];                
    pc = &c[bs];                

    

    for (l=0;l!=n;l++){                                      
      for (j=0;j!=mm;j++) pc[j] = pb[j];
      pc+=ldc;pb+=ldb;                                       
    }                                                        
  }

  

  for (i=mb-1;i!=-1;i--) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    mm = rpntr[i+1]-rpntr[i];   
    pb = &b[bs];                
    pc = &c[bs];                

    for (j=je-1;j!=jb-1;j--) {  
      index = bindx[j];             
      if ( index == i + ind_base ){ 
        continue;                   
      } else {
        pval = &val[indx[j]];       
        ptmp = pval;                                         
                                 
        cs = cpntr[index];          
        cr = cpntr[index+1] - cs;   
        pc = &c[cs];                
        cj = &c[bs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=cr;ii++) {      
            t = 0;
            for (jj=0;jj!=mm;jj++) {        
              t += (*pval++) * cj[jj];      
            }
            pc[ii] -= t;       
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }
  }

  
  

        




}

void VBR_MatTriangSlvUU_CATB_double(
                 const int mb, const int n, 
                  
                 
                 const double *val, const int *indx,
                 const int *bindx, const int *rpntr,
                 const int *cpntr, const int *bpntrb, const int *bpntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc,
                  const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,mm,bs,br,cs,cr;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=rpntr[mb]-rpntr[0];
  double t;

  val-=ind_base;
  b-=ind_base;
  c-=ind_base;
  indx-=ind_base;
  bindx-=ind_base;
  cpntr-=ind_base;



  
  
  mm=0; 
  ds=0;
  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    ds += mm*mm;                
    mm = rpntr[i+1]-rpntr[i];   
    pb = &b[bs];                
    pc = &c[bs];                

    

    for (l=0;l!=n;l++){                                      
      for (j=0;j!=mm;j++) pc[j] = pb[j];
      pc+=ldc;pb+=ldb;                                       
    }                                                        
  }

  

  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    mm = rpntr[i+1]-rpntr[i];   

    for (j=jb;j!=je;j++) {      
      index = bindx[j];             
      if ( index == i + ind_base ){ 
        continue;                   
      } else {
        pval = &val[indx[j]];       
        ptmp = pval;                                         
                                 
        cs = cpntr[index];          
        cr = cpntr[index+1] - cs;   
        pc = &c[cs];                
        cj = &c[bs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=cr;ii++) {      
            t = 0;
            for (jj=0;jj!=mm;jj++) {        
              t += (*pval++) * cj[jj];      
            }
            pc[ii] -= t;
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }
  }

  
  

        




}



void VBR_MatTriangSlvLU_CaAB_double(
                 const int mb, const int n, 
                  
                 const double alpha,
                 const double *val, const int *indx,
                 const int *bindx, const int *rpntr,
                 const int *cpntr, const int *bpntrb, const int *bpntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc,
                  const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,mm,bs,br,cs,cr;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=rpntr[mb]-rpntr[0];
  double t;

  val-=ind_base;
  b-=ind_base;
  c-=ind_base;
  indx-=ind_base;
  bindx-=ind_base;
  cpntr-=ind_base;



  
  
  mm=0; 
  ds=0;
  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    ds += mm*mm;                
    mm = rpntr[i+1]-rpntr[i];   
    pb = &b[bs];                
    pc = &c[bs];                

    

    for (l=0;l!=n;l++){                                      
      for (j=0;j!=mm;j++) pc[j] = pb[j];
      pc+=ldc;pb+=ldb;                                       
    }                                                        
  }

  

  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = rpntr[i];              
    mm = rpntr[i+1]-rpntr[i];   

    for (j=jb;j!=je;j++) {      
      index = bindx[j];             
      if ( index == i + ind_base ){ 
        continue;                   
      } else {
        pval = &val[indx[j]];       
        ptmp = pval;                                         
                                 
        cs = cpntr[index];          
        cr = cpntr[index+1] - cs;   
        pb = &b[bs];                
        pc = &c[bs];                
        cj = &c[cs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=cr;ii++) {      
            t = cj[ii];
            for (jj=0;jj!=mm;jj++) {        
              pc[jj] -= (*pval++) * t;       
            }
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }