dcscvtsl.c

来自「用于求解大型稀疏线性方程组Ax=b的数值计算库.」· C语言 代码 · 共 2,462 行 · 第 1/5 页

}

void CSC_VecTriangSlvUD_CDAB_double(
                 const int m,  const double *dvl,
                   const double *val,
                 const int *indx, const int *pntrb, const int *pntre,
                 const double *b,  
                 double *c,  double *work, const int ind_base)
{
  int i, j, jb, je, index;
  double *pc=c;
  double z; 
  val-=ind_base;
  indx-=ind_base;

    for (i=0;i!=m;i++){
      *pc = b[i];
      pc++;
    }                                     

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

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

    for (i=0;i!=m;i++){
      *pc = b[i];
      pc++;
    }                                     

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




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

    for (i=0;i!=m;i++){
      *pc = b[i];
      pc++;
    }                                     

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

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

    for (i=0;i!=m;i++){
      *pc = b[i];
      pc++;
    }                                     

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

void CSC_VecTriangSlvUD_CaDAB_double(
                 const int m,  const double *dvl,
                  const double alpha, const double *val,
                 const int *indx, const int *pntrb, const int *pntre,
                 const double *b,  
                 double *c,  double *work, const int ind_base)
{
  int i, j, jb, je, index;
  double *pc=c;
  double z; 
  val-=ind_base;
  indx-=ind_base;

    for (i=0;i!=m;i++){
      *pc = b[i];
      pc++;
    }                                     

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

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

    for (i=0;i!=m;i++){
      *pc = b[i];
      pc++;
    }                                     

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




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

    for (i=0;i!=m;i++){
      *pwork++ =  (*pc);                    
      *pc = b[i];
      pc++;
    }                                     

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

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

    for (i=0;i!=m;i++){
      *pwork++ =  (*pc);                    
      *pc = b[i];
      pc++;
    }                                     

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

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

    for (i=0;i!=m;i++){
      *pwork++ =  (*pc);                    
      *pc = b[i];
      pc++;
    }                                     

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

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

    for (i=0;i!=m;i++){
      *pwork++ =  (*pc);                    
      *pc = b[i];
      pc++;
    }                                     

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




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

    for (i=0;i!=m;i++){
      *pwork++ =  (*pc);                    
      *pc = b[i];
      pc++;
    }                                     

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

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

    for (i=0;i!=m;i++){
      *pwork++ =  (*pc);                    
      *pc = b[i];
      pc++;
    }                                     

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

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

    for (i=0;i!=m;i++){
      *pwork++ =  (*pc);                    
      *pc = b[i];
      pc++;
    }                                     

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

void CSC_VecTriangSlvUU_CaDABC_double(