syminitvectors.c

数学算法的实现库。可以实现常见的线性计算。

    // remember that scaling is explicitly applied to A
    scale=PRE.colscal;


    // seek for an exact solution
    if (rhstyp[2]=='X' || rhstyp[2]=='x') {
       j=nrhs*n;
       if (rhstyp[1]=='G' || rhstyp[1]=='g') 
	 j*=2;
       if (rhstyp[0]=='M' || rhstyp[0]=='m')
	 j-=nrhs*n;
       for (i=0; i<n; i++,j++)
	   sol[i+A.nr*l]=rhs[j+A.nr*l];
    }


    // release part of rhs that may store the uncompressed rhs
    if (nrhs!=0 && (rhstyp[0]=='M' || rhstyp[0]=='m')) {
       rhs-=n;
       rhs-=A.nr*l;
    }

    // set up right hand side
    if (nrhs==0) {
      // prescribe artificial solution
      for (i=0; i<n; i++)  {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	  sol[i+A.nr*l]=1.0+i*l;
#else
	  sol[i+A.nr*l].r=1.0+i*(l-1);
	  sol[i+A.nr*l].i=0.1*(A.nr-i+l-1);
#endif
      } // end for i

      // construct associated right hand side rhs=A*sol
      // remember that A is already rescaled!!!
      // rescale solution
      for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_ 
	  sol[i+A.nr*l]/=scale[i];
#else
	  val=1.0/(scale[i].r*scale[i].r +
		   scale[i].i*scale[i].i);
	  vb=sol[i+A.nr*l].r;
	  sol[i+A.nr*l].r=( vb*scale[i].r+sol[i+A.nr*l].i*scale[i].i)*val;
	  sol[i+A.nr*l].i=(-vb*scale[i].i+sol[i+A.nr*l].i*scale[i].r)*val;
#endif
      } // end for i

      MYSMATVEC(A,sol+A.nr*l,rhs+A.nr*l);

      // rescale right hand side and solution back to their original
      // counterparts associated with the input matrix
      for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_ 
          sol[i+A.nr*l]*=scale[i];
	  rhs[i+A.nr*l]/=scale[i];
#else
	  vb=sol[i+A.nr*l].r;
	  sol[i+A.nr*l].r=vb*scale[i].r-sol[i+A.nr*l].i*scale[i].i;
	  sol[i+A.nr*l].i=vb*scale[i].i+sol[i+A.nr*l].i*scale[i].r;
	  val=1.0/(scale[i].r*scale[i].r +
		   scale[i].i*scale[i].i);
	  vb=rhs[i+A.nr*l].r;
#ifdef _COMPLEX_SYMMETRIC_
	  // transposed scaling
	  rhs[i+A.nr*l].r=( vb*scale[i].r+rhs[i+A.nr*l].i*scale[i].i)*val;
	  rhs[i+A.nr*l].i=(-vb*scale[i].i+rhs[i+A.nr*l].i*scale[i].r)*val;
#else
	  // conjugate transposed scaling
	  rhs[i+A.nr*l].r=( vb*scale[i].r-rhs[i+A.nr*l].i*scale[i].i)*val;
	  rhs[i+A.nr*l].i=( vb*scale[i].i+rhs[i+A.nr*l].i*scale[i].r)*val;
#endif
#endif
      } // end for i
    } // end if (nrhs==0)
    else {
       if (rhstyp[0]=='M' || rhstyp[0]=='m') {
	  for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_ 
	      rhs[i+A.nr*l]=0;
#else
	      rhs[i+A.nr*l].r=rhs[i+A.nr*l].i=0;
#endif
	  }// end for i

	  // uncompress rhs
	  for (i=0; i<rhsptr[l+1]-rhsptr[l]; i++) {
	      j=rhsind[i]-1;
	      rhs[j+A.nr*l]=rhsval[i];
	  } // end for i
       } // end if
    } // end if-else


    // initial solution
    if (rhstyp[1]=='G' || rhstyp[1]=='g') {
       j=nrhs*n;
       if (rhstyp[0]=='M' || rhstyp[0]=='m')
	  j=n;
       for (i=0; i<n; i++,j++)
	   sol[i+A.nr*l]=rhs[j+A.nr*l];
    }
    else
       for (i=0; i<n; i++)
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	   sol[i+A.nr*l]=0;
#else
           sol[i+A.nr*l].r=sol[i+A.nr*l].i=0;
#endif

    // ------   compute initial residual   -----
    // remember that A is already rescaled!!!

    // rescale solution
    for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_ 
        sol[i+A.nr*l]/=scale[i];
#else
	val=1.0/(scale[i].r*scale[i].r +
  	         scale[i].i*scale[i].i);
	vb=sol[i+A.nr*l].r;
	sol[i+A.nr*l].r=( vb*scale[i].r+sol[i+A.nr*l].i*scale[i].i)*val;
	sol[i+A.nr*l].i=(-vb*scale[i].i+sol[i+A.nr*l].i*scale[i].r)*val;
#endif
    } // end for i

    MYSMATVEC(A,sol+A.nr*l,param.dbuff);

    for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_ 
        sol[i+A.nr*l]*=scale[i];
#else
	vb=sol[i+A.nr*l].r;
	sol[i+A.nr*l].r=vb*scale[i].r-sol[i+A.nr*l].i*scale[i].i;
	sol[i+A.nr*l].i=vb*scale[i].i+sol[i+A.nr*l].i*scale[i].r;
#endif
    } // end for i

    for (i=0; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_ 
        param.dbuff[i]-=rhs[i+A.nr*l]*scale[i];
#else
#ifdef _COMPLEX_SYMMETRIC_
	// transposed scaling
	param.dbuff[i].r-=rhs[i+A.nr*l].r*scale[i].r-rhs[i+A.nr*l].i*scale[i].i;
	param.dbuff[i].i-=rhs[i+A.nr*l].r*scale[i].i+rhs[i+A.nr*l].i*scale[i].r;
#else
	// conjugate transposed scaling
	param.dbuff[i].r-= rhs[i+A.nr*l].r*scale[i].r+rhs[i+A.nr*l].i*scale[i].i;
	param.dbuff[i].i-=-rhs[i+A.nr*l].r*scale[i].i+rhs[i+A.nr*l].i*scale[i].r;
#endif
#endif
    } // end for i

    i=1;
    val=NRM(&n, param.dbuff, &i);
    // -----------------------------------------