dbsrvtsl.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:35:50 EDT 1996 */

#include "dbsrvtsl.h"



void BSR_VecTriangSlvLU_CAB_double(
                 const int mb,  
                  
                 
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b,   
                 double *c, 
                  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,bs,br,cs;
  int index;
  int ind;
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;



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

    

      for (j=0;j!=lb;j++) pc[j] = pb[j];
  }

  

  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = i*lb;                  

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

  
  

        



}


void BSR_VecTriangSlvUU_CAB_double(
                 const int mb, 
                  
                 
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b,   
                 double *c, 
                  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,bs,br,cs;
  int index;
  int ind;
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;



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

    

      for (j=0;j!=lb;j++) pc[j] = pb[j];
  }

  

  for (i=mb-1;i!=-1;i--) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = lb*i;                  

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

  
  

        




}

void BSR_VecTriangSlvLU_CATB_double(
                 const int mb, 
                  
                 
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b,   
                 double *c, 
                  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,bs,br,cs;
  int index;
  int ind;
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;



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

    

      for (j=0;j!=lb;j++) pc[j] = pb[j];
  }

  

  for (i=mb-1;i!=-1;i--) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = i*lb;                  

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

  
  

        




}

void BSR_VecTriangSlvUU_CATB_double(
                 const int mb,  
                  
                 
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b,   
                 double *c, 
                  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,bs,br,cs;
  int index;
  int ind;
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;



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

    

      for (j=0;j!=lb;j++) pc[j] = pb[j];
  }

  

  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = i*lb;                  

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

  
  

        




}



void BSR_VecTriangSlvLU_CaAB_double(
                 const int mb,  
                  
                 const double alpha,
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b,   
                 double *c, 
                  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,bs,br,cs;
  int index;
  int ind;
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;



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

    

      for (j=0;j!=lb;j++) pc[j] = pb[j];
  }

  

  for (i=0;i!=mb;i++) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = i*lb;                  

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

  
  

        

  
  

  pc=c;                                   
    for (i=0;i!=m;i++) {                  
          *pc *= alpha;                   
          pc++;                           
    }                                     

}


void BSR_VecTriangSlvUU_CaAB_double(
                 const int mb, 
                  
                 const double alpha,
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b,   
                 double *c, 
                  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,bs,br,cs;
  int index;
  int ind;
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;



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

    

      for (j=0;j!=lb;j++) pc[j] = pb[j];
  }

  

  for (i=mb-1;i!=-1;i--) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = lb*i;                  

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

  
  

        

  
  
  

  pc=c;                                   
    for (i=0;i!=m;i++) {                  
          *pc *= alpha;                   
          pc++;                           
    }                                     


}

void BSR_VecTriangSlvLU_CaATB_double(
                 const int mb, 
                  
                 const double alpha,
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b,   
                 double *c, 
                  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,bs,br,cs;
  int index;
  int ind;
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;



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

    

      for (j=0;j!=lb;j++) pc[j] = pb[j];
  }

  

  for (i=mb-1;i!=-1;i--) {      
    jb = bpntrb[i];             
    je = bpntre[i];             
    bs = i*lb;                  

    for (j=je-1;j!=jb-1;j--) {  
      index = bindx[j];             
      if ( index == i + ind_base ){ 
        continue;                   
      } else {
        pval = &val[(j-1)*mm];      
                                 
        cs = (index-1)*lb;
        pc = &c[cs];                
        cj = &c[bs];