dcsrmtsl.c

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

    }
    for (i=0;i!=m;i++) {                  
        *pc *= alpha * dvl[i];            
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}




void CSR_MatTriangSlvLD_CDABC_double(
                 const int m, const int n, const double *dvl, 
                   const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je, index;
  double *pc=c;
  double *pwork=work;                             
  int l;                                              
  double z; 
  double valtmp;
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=0;i!=m;i++) {
      z = 0;
      jb = pntrb[i];
      je = pntre[i];
      for (j=jb;j<je;j++) {
        index = indx[j];
        if ( index == i+ind_base ) {
          valtmp = val[j];
        } else {
          z += c[index] * val[j];
        }
      }
      pc[i] = (b[i] - z) / valtmp;
    }
    for (i=0;i!=m;i++) {                  
        *pc *=  dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}

void CSR_MatTriangSlvLU_CDABC_double(
                 const int m, const int n, const double *dvl, 
                   const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je;
  double *pc=c;
  double *pwork=work;                                
  int l;                                              
  double z; 
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=0;i!=m;i++) {
      z = 0;
      jb = pntrb[i];
      je = pntre[i];
      for (j=jb;j<je;j++) {
        z += c[indx[j]] * val[j];
      }
      pc[i] = b[i] - z;
    }
    for (i=0;i!=m;i++) {                  
        *pc *=  dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}

void CSR_MatTriangSlvUD_CDABC_double(
                 const int m, const int n, const double *dvl, 
                   const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je, index;
  double *pc=c;
  double *pwork=work;                                
  double valtmp;
  int l;                                              
  double z; 
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=m-1;i!=-1; i--) {
      z = 0;
      jb = pntrb[i];
      je =  pntre[i];
      for (j=jb+1; j<je; j++) {
        z += c[indx[j]] * val[j];
      }
      pc[i] = (b[i] - z) / val[jb];
    }
    for (i=0;i!=m;i++) {                  
        *pc *=  dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}

void CSR_MatTriangSlvUU_CDABC_double(
                 const int m, const int n, const double *dvl, 
                   const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je;
  double *pc=c;
  double *pwork=work;                                
  int l;                                              
  double z; 
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=m-1;i!=-1;i--) {
      z = 0;
      jb =  pntrb[i];
      je =  pntre[i];
      for (j=jb; j<je;j++) {
        z += c[indx[j]] * val[j];
      }
      pc[i] = b[i] - z;
    }
    for (i=0;i!=m;i++) {                  
        *pc *=  dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}




void CSR_MatTriangSlvLD_CaDABC_double(
                 const int m, const int n, const double *dvl, 
                  const double alpha, const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je, index;
  double *pc=c;
  double *pwork=work;                             
  int l;                                              
  double z; 
  double valtmp;
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=0;i!=m;i++) {
      z = 0;
      jb = pntrb[i];
      je = pntre[i];
      for (j=jb;j<je;j++) {
        index = indx[j];
        if ( index == i+ind_base ) {
          valtmp = val[j];
        } else {
          z += c[index] * val[j];
        }
      }
      pc[i] = (b[i] - z) / valtmp;
    }
    for (i=0;i!=m;i++) {                  
        *pc *= alpha * dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}

void CSR_MatTriangSlvLU_CaDABC_double(
                 const int m, const int n, const double *dvl, 
                  const double alpha, const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je;
  double *pc=c;
  double *pwork=work;                                
  int l;                                              
  double z; 
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=0;i!=m;i++) {
      z = 0;
      jb = pntrb[i];
      je = pntre[i];
      for (j=jb;j<je;j++) {
        z += c[indx[j]] * val[j];
      }
      pc[i] = b[i] - z;
    }
    for (i=0;i!=m;i++) {                  
        *pc *= alpha * dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}

void CSR_MatTriangSlvUD_CaDABC_double(
                 const int m, const int n, const double *dvl, 
                  const double alpha, const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je, index;
  double *pc=c;
  double *pwork=work;                                
  double valtmp;
  int l;                                              
  double z; 
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=m-1;i!=-1; i--) {
      z = 0;
      jb = pntrb[i];
      je =  pntre[i];
      for (j=jb+1; j<je; j++) {
        z += c[indx[j]] * val[j];
      }
      pc[i] = (b[i] - z) / val[jb];
    }
    for (i=0;i!=m;i++) {                  
        *pc *= alpha * dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}

void CSR_MatTriangSlvUU_CaDABC_double(
                 const int m, const int n, const double *dvl, 
                  const double alpha, const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb,  
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je;
  double *pc=c;
  double *pwork=work;                                
  int l;                                              
  double z; 
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=m-1;i!=-1;i--) {
      z = 0;
      jb =  pntrb[i];
      je =  pntre[i];
      for (j=jb; j<je;j++) {
        z += c[indx[j]] * val[j];
      }
      pc[i] = b[i] - z;
    }
    for (i=0;i!=m;i++) {                  
        *pc *= alpha * dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}




void CSR_MatTriangSlvLD_CDABbC_double(
                 const int m, const int n, const double *dvl, 
                   const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb, const double beta, 
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je, index;
  double *pc=c;
  double *pwork=work;                             
  int l;                                              
  double z; 
  double valtmp;
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {                                
    pc=c+ind_base;
    for (i=0;i!=m;i++) {
      z = 0;
      jb = pntrb[i];
      je = pntre[i];
      for (j=jb;j<je;j++) {
        index = indx[j];
        if ( index == i+ind_base ) {
          valtmp = val[j];
        } else {
          z += c[index] * val[j];
        }
      }
      pc[i] = (b[i] - z) / valtmp;
    }
    for (i=0;i!=m;i++) {                  
        *pc *=  dvl[i];            
        *pc += *pwork;                    
        pwork++;                          
        pc++;                             
    }                                     
    c += ldc; b += ldb;                               
  }                                                   
}

void CSR_MatTriangSlvLU_CDABbC_double(
                 const int m, const int n, const double *dvl, 
                   const double *val, 
                 const int *indx, const int *pntrb, const int *pntre, 
                 const double *b, const int ldb, const double beta, 
                 double *c, const int ldc, double *work, const int ind_base)
{
  int i, j, jb, je;
  double *pc=c;
  double *pwork=work;                                
  int l;                                              
  double z; 
  val-=ind_base;
  indx-=ind_base;

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

  pwork=work;                                        
  c-=ind_base;
  for (l=0;l!=n;l++) {