dbsrmtsl.c

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

  const double *pdr;
  const double *pdl;
  double *pwork=work;                             
  double *ptmpwork;                               
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,bs,br,cs;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;

  
  
  

  for (l=0;l!=n;l++)                              
    for (i=0;i!=m;i++) *pwork++ =  (*pc++);       

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

    

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

  

  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];      
        ptmp = pval;                                         
                                 
        cs = (index-1)*lb;          
        pc = &c[cs];                
        cj = &c[bs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=lb;ii++) {      
            t = 0;
            for (jj=0;jj!=lb;jj++) {        
              t += (*pval++) * cj[jj];      
            }
            pc[ii] -= t;
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }
  }

  
  

        

  
  
  

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


}



void BSR_MatTriangSlvLU_CABbC_double(
                 const int mb, const int n, 
                  
                 
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b, const int ldb, const double beta, 
                 double *c, const int ldc,
                 double *work, const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  double *pwork=work;                             
  double *ptmpwork;                               
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,bs,br,cs;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;

  
  
  

  for (l=0;l!=n;l++)                              
    for (i=0;i!=m;i++) *pwork++ = beta * (*pc++);       

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

    

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

  

  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];      
        ptmp = pval;                                         
                                 
        cs = (index-1)*lb;          
        pc = &c[bs];                
        cj = &c[cs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=lb;ii++) {      
            t = cj[ii];
            for (jj=0;jj!=lb;jj++) {        
              pc[jj] -= (*pval++) * t;       
            }
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }
  }

  
  

        

  
  

  pc=c;                                   
  pwork=work;                             
  for (l=0;l!=n;l++)                      
    for (i=0;i!=m;i++) {                  
          *pc += *pwork;                  
          pc++;                           
          pwork++;                        
    }                                     

}


void BSR_MatTriangSlvUU_CABbC_double(
                 const int mb, const int n,
                  
                 
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b, const int ldb, const double beta, 
                 double *c, const int ldc,
                 double *work, const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  double *pwork=work;                            
  double *ptmpwork;                              
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,bs,br,cs;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;

  
  
  

  for (l=0;l!=n;l++)                              
    for (i=0;i!=m;i++) *pwork++ = beta * (*pc++);       

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

    

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

  

  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];      
        ptmp = pval;                                         
                                 
        cs = (index-1)*lb;          
        pc = &c[bs];                
        cj = &c[cs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=lb;ii++) {      
            t = cj[ii];
            for (jj=0;jj!=lb;jj++) {        
              pc[jj] -= (*pval++) * t;       
            }
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }
  }

  
  

        

  
  
  

  pc=c;                                   
  pwork=work;                             
  for (l=0;l!=n;l++)                      
    for (i=0;i!=m;i++) {                  
          *pc += *pwork;                  
          pc++;                           
          pwork++;                        
    }                                     


}

void BSR_MatTriangSlvLU_CATBbC_double(
                 const int mb, const int n,
                  
                 
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b, const int ldb, const double beta, 
                 double *c, const int ldc,
                 double *work, const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  double *pwork=work;                             
  double *ptmpwork;                               
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,bs,br,cs;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;

  
  
  

  for (l=0;l!=n;l++)                              
    for (i=0;i!=m;i++) *pwork++ = beta * (*pc++);       

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

    

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

  

  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];      
        ptmp = pval;                                         
                                 
        cs = (index-1)*lb;
        pc = &c[cs];                
        cj = &c[bs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=lb;ii++) {      
            t = 0;
            for (jj=0;jj!=lb;jj++) {        
              t += (*pval++) * cj[jj];      
            }
            pc[ii] -= t;       
          }
          pc+=ldc;cj+=ldc;pval=ptmp; 
        }                                                    
      }
    }
  }

  
  

        

  
  
  

  pc=c;                                   
  pwork=work;                             
  for (l=0;l!=n;l++)                      
    for (i=0;i!=m;i++) {                  
          *pc += *pwork;                  
          pc++;                           
          pwork++;                        
    }                                     


}

void BSR_MatTriangSlvUU_CATBbC_double(
                 const int mb, const int n, 
                  
                 
                 const double *val, const int *bindx, 
                 const int *bpntrb, const int *bpntre, const int lb,
                 const double *b, const int ldb, const double beta, 
                 double *c, const int ldc,
                 double *work, const int ind_base)
{
  const double *pb;
  double *pc=c;
  double *cj;
  const double *pdr;
  const double *pdl;
  double *pwork=work;                             
  double *ptmpwork;                               
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int ds,bs,br,cs;
  int index;
  int ind;
  int l;                                               
  int ii,jj;
  int m=mb*lb;
  int mm=lb*lb;
  double t;

  bindx-=ind_base;

  
  
  

  for (l=0;l!=n;l++)                              
    for (i=0;i!=m;i++) *pwork++ = beta * (*pc++);       

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

    

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

  

  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];      
        ptmp = pval;                                         
                                 
        cs = (index-1)*lb;          
        pc = &c[cs];                
        cj = &c[bs];                
        for (l=0;l!=n;l++) {            
          for (ii=0;ii!=lb;ii++) {      
            t = 0;
            for (jj=0;jj!=lb;jj++) {        
              t += (*pval++) * cj[jj];      
            }
            pc[ii] -= t;
          }