dvbrmtsl.c

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

  ptmpwork = pwork;          
                             
                             
                             
                             
                             
  
  
  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; 
        }                                                    
      }
    }
  }

  
  

        

  
  
  

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


}

void VBR_MatTriangSlvLU_CATBC_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,
                 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,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;

  
  
  

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

  ptmpwork = pwork;        
                           
                           
                           
                           
                           
  
  
  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; 
        }                                                    
      }
    }
  }

  
  

        

  
  
  

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


}

void VBR_MatTriangSlvUU_CATBC_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,
                 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,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;

  
  
  

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

  ptmpwork = pwork;          
                             
                             
                             
                             
                             
  
  
  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; 
        }                                                    
      }
    }
  }

  
  

        

  
  
  

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


}



void VBR_MatTriangSlvLU_CaABC_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,
                 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,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;

  
  
  

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

  ptmpwork = pwork;          
                             
                             
                             
                             
                             
  
  
  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; 
        }                                                    
      }
    }
  }

  
  

        

  
  
  

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

}


void VBR_MatTriangSlvUU_CaABC_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,
                 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,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;

  
  
  

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

  ptmpwork = pwork;          
                             
                             
                             
                             
                             
  
  
  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];