dbcovml.c

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

  }
}




void BCO_VecMult_CABbC_double(
                 const int mb,  const int kb, 
                 const double *val, const int *bindx, const int *bjndx,
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  int i,j,jb,je,block;
  int cs,bs,br;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;

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

  for (i=0;i!=bnnz;i++) {      
      cs = (bindx[i]-1)*lb;     
      bs = (bjndx[i]-1)*lb;        
      pval = &val[i*mm];           
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
         if( pb[jj] != 0.0 ) {           
           t =  pb[jj];            
           for (ii=0;ii!=lb;ii++) {      
             pc[ii] += t* (*pval++);        
           }
         } else {
           pval+=lb;                     
         }
       }
  }
}

void BCO_VecMult_CATBbC_double(
                 const int mb,  const int kb, 
                 const double *val, const int *bindx, const int *bjndx, 
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  int i,j,jb,je,block;
  int cs,bs,cr;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;

    for (i=0;i!=k;i++) *pc++ *= beta;                  

  for (i=0;i!=bnnz;i++) {      
      cs = (bjndx[i]-1)*lb;     
      bs = (bindx[i]-1)*lb;        
      pval = &val[i*mm];           
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
           t = 0;
           for (ii=0;ii!=lb;ii++) {      
             t +=  pb[ii] * (*pval++);
           }
           pc[jj] += t;                  
       }
  } }

void BCOsymm_VecMult_CABbC_double(
                 const int mb,  const int kb, 
                 const double *val, const int *bindx, const int *bjndx, 
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  const double *ptmp; 
  int i,j,jb,je,block;
  int cs,css,cr,bs,br;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;

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

  for (i=0;i!=bnnz;i++) {    
      cs = (bindx[i]-1)*lb;     
      bs = (bjndx[i]-1)*lb;        
      pval = &val[i*mm];           
      ptmp = pval;                 
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
         if( pb[jj] != 0.0 ) {           
           t =  pb[jj];            
           for (ii=0;ii!=lb;ii++) {      
             pc[ii] += t* (*pval++);        
           }
         } else {
           pval+=lb;                     
         }
       }
 
      if ( cs == bs ) {
        continue;                  
      }
 
      pval = ptmp;                 
                                   
      pb = &b[cs];                 
      pc = &c[bs];                 
       for (jj=0;jj!=lb;jj++) {       
           t = 0;
           for (ii=0;ii!=lb;ii++) {      
             t +=  pb[ii] * (*pval++);
           }
           pc[jj] += t;                  
       }
  }
}

void BCOskew_VecMult_CABbC_double(
                 const int mb,  const int kb, 
                 const double *val, const int *bindx, const int *bjndx, 
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int cs,css,cr,bs,br;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;

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

  for (i=0;i!=bnnz;i++) {    
      cs = (bindx[i]-1)*lb;     
      bs = (bjndx[i]-1)*lb;        
 
      if ( cs == bs ) {
        continue;                  
      }
      pval = &val[i*mm];           
      ptmp = pval;                 
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
         if( pb[jj] != 0.0 ) {           
           t =  pb[jj];            
           for (ii=0;ii!=lb;ii++) {      
             pc[ii] += t* (*pval++);        
           }
         } else {
           pval+=lb;                     
         }
       }
 
      pval = ptmp;                 
                                   
      pb = &b[cs];                 
      pc = &c[bs];                 
       for (jj=0;jj!=lb;jj++) {       
           t = 0;
           for (ii=0;ii!=lb;ii++) {      
             t -=  pb[ii] * (*pval++);
           }
           pc[jj] += t;                  
       }
  }
}

void BCOskew_VecMult_CATBbC_double(
                 const int mb,  const int kb, 
                 const double *val, const int *bindx, const int *bjndx,
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  int i,j,jb,je,block;
  int cs,css,cr,bs,br;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;
  const double *ptmp;

    for (i=0;i!=k;i++) *pc++ *= beta;                  

  for (i=0;i!=bnnz;i++) {    
      cs = (bindx[i]-1)*lb;     
      bs = (bjndx[i]-1)*lb;        
 
      if ( cs == bs ) {
        continue;                  
      }
      pval = &val[i*mm];           
      ptmp = pval;                 
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
         if( pb[jj] != 0.0 ) {           
           t =  pb[jj];            
           for (ii=0;ii!=lb;ii++) {      
             pc[ii] -= t* (*pval++);        
           }
         } else {
           pval+=lb;                     
         }
       }
 
      pval = ptmp;                 
                                   
      pb = &b[cs];                 
      pc = &c[bs];                 
       for (jj=0;jj!=lb;jj++) {       
           t = 0;
           for (ii=0;ii!=lb;ii++) {      
             t +=  pb[ii] * (*pval++);
           }
           pc[jj] += t;                  
       }
  }
}




void BCO_VecMult_CaABbC_double(
                 const int mb,  const int kb, const double alpha,
                 const double *val, const int *bindx, const int *bjndx,
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  int i,j,jb,je,block;
  int cs,bs,br;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;

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

  for (i=0;i!=bnnz;i++) {      
      cs = (bindx[i]-1)*lb;     
      bs = (bjndx[i]-1)*lb;        
      pval = &val[i*mm];           
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
         if( pb[jj] != 0.0 ) {           
           t = alpha * pb[jj];            
           for (ii=0;ii!=lb;ii++) {      
             pc[ii] += t* (*pval++);        
           }
         } else {
           pval+=lb;                     
         }
       }
  }
}

void BCO_VecMult_CaATBbC_double(
                 const int mb,  const int kb, const double alpha,
                 const double *val, const int *bindx, const int *bjndx, 
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  int i,j,jb,je,block;
  int cs,bs,cr;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;

    for (i=0;i!=k;i++) *pc++ *= beta;                  

  for (i=0;i!=bnnz;i++) {      
      cs = (bjndx[i]-1)*lb;     
      bs = (bindx[i]-1)*lb;        
      pval = &val[i*mm];           
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
           t = 0;
           for (ii=0;ii!=lb;ii++) {      
             t += alpha * pb[ii] * (*pval++);
           }
           pc[jj] += t;                  
       }
  } }

void BCOsymm_VecMult_CaABbC_double(
                 const int mb,  const int kb, const double alpha,
                 const double *val, const int *bindx, const int *bjndx, 
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  const double *ptmp; 
  int i,j,jb,je,block;
  int cs,css,cr,bs,br;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;

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

  for (i=0;i!=bnnz;i++) {    
      cs = (bindx[i]-1)*lb;     
      bs = (bjndx[i]-1)*lb;        
      pval = &val[i*mm];           
      ptmp = pval;                 
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
         if( pb[jj] != 0.0 ) {           
           t = alpha * pb[jj];            
           for (ii=0;ii!=lb;ii++) {      
             pc[ii] += t* (*pval++);        
           }
         } else {
           pval+=lb;                     
         }
       }
 
      if ( cs == bs ) {
        continue;                  
      }
 
      pval = ptmp;                 
                                   
      pb = &b[cs];                 
      pc = &c[bs];                 
       for (jj=0;jj!=lb;jj++) {       
           t = 0;
           for (ii=0;ii!=lb;ii++) {      
             t += alpha * pb[ii] * (*pval++);
           }
           pc[jj] += t;                  
       }
  }
}

void BCOskew_VecMult_CaABbC_double(
                 const int mb,  const int kb, const double alpha,
                 const double *val, const int *bindx, const int *bjndx, 
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  const double *ptmp;
  int i,j,jb,je,block;
  int cs,css,cr,bs,br;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;

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

  for (i=0;i!=bnnz;i++) {    
      cs = (bindx[i]-1)*lb;     
      bs = (bjndx[i]-1)*lb;        
 
      if ( cs == bs ) {
        continue;                  
      }
      pval = &val[i*mm];           
      ptmp = pval;                 
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
         if( pb[jj] != 0.0 ) {           
           t = alpha * pb[jj];            
           for (ii=0;ii!=lb;ii++) {      
             pc[ii] += t* (*pval++);        
           }
         } else {
           pval+=lb;                     
         }
       }
 
      pval = ptmp;                 
                                   
      pb = &b[cs];                 
      pc = &c[bs];                 
       for (jj=0;jj!=lb;jj++) {       
           t = 0;
           for (ii=0;ii!=lb;ii++) {      
             t -= alpha * pb[ii] * (*pval++);
           }
           pc[jj] += t;                  
       }
  }
}

void BCOskew_VecMult_CaATBbC_double(
                 const int mb,  const int kb, const double alpha,
                 const double *val, const int *bindx, const int *bjndx,
                 const int bnnz, const int lb,
                 const double *b,  const double beta, 
                 double *c, 
                 const int ind_base)
{
  double t;
  const double *pb;
  double *pc=c;
  const double *pval;
  int i,j,jb,je,block;
  int cs,css,cr,bs,br;
  int index;
  int ii,jj;
  int m=mb*lb;
  int k=kb*lb;
  int mm=lb*lb;
  const double *ptmp;

    for (i=0;i!=k;i++) *pc++ *= beta;                  

  for (i=0;i!=bnnz;i++) {    
      cs = (bindx[i]-1)*lb;     
      bs = (bjndx[i]-1)*lb;        
 
      if ( cs == bs ) {
        continue;                  
      }
      pval = &val[i*mm];           
      ptmp = pval;                 
                                   
      pb = &b[bs];                 
      pc = &c[cs];                 
       for (jj=0;jj!=lb;jj++) {       
         if( pb[jj] != 0.0 ) {           
           t = alpha * pb[jj];            
           for (ii=0;ii!=lb;ii++) {      
             pc[ii] -= t* (*pval++);        
           }
         } else {
           pval+=lb;                     
         }
       }
 
      pval = ptmp;                 
                                   
      pb = &b[cs];                 
      pc = &c[bs];                 
       for (jj=0;jj!=lb;jj++) {       
           t = 0;
           for (ii=0;ii!=lb;ii++) {      
             t += alpha * pb[ii] * (*pval++);
           }
           pc[jj] += t;                  
       }
  }
}