dvbrvtsl.c

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

                 const int mb, 
                  
                 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 double beta, 
                 double *c, 
                 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 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 (i=0;i!=m;i++) *pwork++ = beta * (*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 (j=0;j!=mm;j++) pc[j] = pb[j];
  }

  

  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]];       
                                 
        cs = cpntr[index];          
        cr = cpntr[index+1] - cs;   
        pb = &b[bs];                
        pc = &c[bs];                
        cj = &c[cs];                
          for (ii=0;ii!=cr;ii++) {      
            t = cj[ii];
            for (jj=0;jj!=mm;jj++) {        
              pc[jj] -= (*pval++) * t;       
            }
          }
      }
    }
  }

  
  

        

  
  
  

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


}

void VBR_VecTriangSlvLU_CaATBbC_double(
                 const int mb, 
                  
                 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 double beta, 
                 double *c, 
                 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 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 (i=0;i!=m;i++) *pwork++ = beta * (*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 (j=0;j!=mm;j++) pc[j] = pb[j];
  }

  

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

  
  

        

  
  
  

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


}

void VBR_VecTriangSlvUU_CaATBbC_double(
                 const int mb,  
                  
                 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 double beta, 
                 double *c, 
                 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 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 (i=0;i!=m;i++) *pwork++ = beta * (*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 (j=0;j!=mm;j++) pc[j] = pb[j];
  }

  

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

  
  

        

  
  
  

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


}



void VBR_VecTriangSlvLU_CDAB_double(
                 const int mb,  
                 const double *dvl, 
                 
                 const double *val, const int *indx,
                 const int *bindx, const int *rpntr,
                 const int *cpntr, const int *bpntrb, const int *bpntre, 
                 const double *b,   
                 double *c, 
                 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 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;



  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 (j=0;j!=mm;j++) pc[j] = pb[j];
  }

  

  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]];       
                                 
        cs = cpntr[index];          
        cr = cpntr[index+1] - cs;   
        pb = &b[bs];                
        pc = &c[bs];                
        cj = &c[cs];                
          for (ii=0;ii!=cr;ii++) {      
            t = cj[ii];
            for (jj=0;jj!=mm;jj++) {        
              pc[jj] -= (*pval++) * t;       
            }
          }
      }
    }
  }

  
  

  ds = 0;                                                    
  mm = 0;                                                    
  for (i=0;i!=mb;i++) {      
    bs = rpntr[i];              
    ds += mm*mm;                
    mm = rpntr[i+1]-rpntr[i];   
    pdl = &dvl[ds];                                          
    pc = &c[bs];                                             
    pwork=ptmpwork;                                          
      for (j=0;j!=mm;j++) {     
        t = 0;                                               
        for (jj=0;jj!=mm;jj++)                               
          t+= pc[jj]*pdl[j+jj*mm];  
        pwork[j] = t;               
      }                                                      
        
      for (j=0;j!=mm;j++) pc[j] = pwork[j];                   
  }                                                          



}


void VBR_VecTriangSlvUU_CDAB_double(
                 const int mb, 
                 const double *dvl, 
                 
                 const double *val, const int *indx,
                 const int *bindx, const int *rpntr,
                 const int *cpntr, const int *bpntrb, const int *bpntre, 
                 const double *b,   
                 double *c, 
                 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 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;



  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 (j=0;j!=mm;j++) pc[j] = pb[j];
  }

  

  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]];       
                                 
        cs = cpntr[index];          
        cr = cpntr[index+1] - cs;   
        pb = &b[bs];                
        pc = &c[bs];                
        cj = &c[cs];                
          for (ii=0;ii!=cr;ii++) {      
            t = cj[ii];
            for (jj=0;jj!=mm;jj++) {        
              pc[jj] -= (*pval++) * t;       
            }
          }
      }
    }
  }

  
  

  ds = 0;                                                    
  mm = 0;                                                    
  for (i=0;i!=mb;i++) {      
    bs = rpntr[i];              
    ds += mm*mm;                
    mm = rpntr[i+1]-rpntr[i];   
    pdl = &dvl[ds];                                          
    pc = &c[bs];                                             
    pwork=ptmpwork;                                          
      for (j=0;j!=mm;j++) {     
        t = 0;                                               
        for (jj=0;jj!=mm;jj++)