unsym.cxx

不错的国外的有限元程序代码,附带详细的manual,可以节省很多的底层工作.

FloatArray*  SkylineUnsym :: backSubstitutionWith (FloatArray* y)
   // Returns the solution x of the system U.x = y , where U is the upper
   // half of the receiver. Note : x overwrites y.
{
   int        i,k,start ;
   double     yK ;
   RowColumn* rowColumnK ;

   for (k=size ; k>0 ; k--) {
      rowColumnK = this->giveRowColumn(k) ;
      yK       = y->at(k) ;
      start    = rowColumnK->giveStart() ;
      for (i=start ; i<k ; i++)
	 y->at(i) -= rowColumnK->atU(i) * yK ;}

   return y ;
}


Skyline*  SkylineUnsym :: diagonalScalingWith (FloatArray* y)
   // Scales y by the diagonal of the receiver. Returns the receiver.
{
   double       diag ;
   int          k ;
   const double precision = 1.e-10 ;

   for (k=1 ; k<=size ; k++) {
      diag = this->giveRowColumn(k)->atDiag() ;
#     ifdef DEBUG
	 if (fabs(diag) < precision) {
	   printf ("error in SkylineUnsym::diagScaling : pivot %d is zero \n",k);
	   exit(0) ;}
#     endif
      y->at(k) /= diag ;}

   return this ;
}


Skyline*  SkylineUnsym :: factorized ()
   // Returns the receiver in L.D.U factorized form. From Golub & van Loan,
   // 1rst edition, pp 83-84.
{
   RowColumn     *rowColumnK,*rowColumnI ;
   FloatArray    *r,*w ;
   double        diag ;
   int           k,i,p,start,startK,startI ;
   const double  PRECISION = 1.e-10 ;

   if (isFactorized)
      return this ;

   if (! size) {
      printf ("error : cannot factorize null-sized matrix\n") ;
      exit(0) ;}

   for (k=1 ; k<=size ; k++) {

      rowColumnK = this -> giveRowColumn(k) ;
      startK   = rowColumnK -> giveStart() ;

      // compute vectors r and w
      r = new FloatArray(k-1) ;
      w = new FloatArray(k-1) ;
      for (p=startK ; p<k ; p++) {
	 diag     = this->giveRowColumn(p)->atDiag() ;
	 r->at(p) = diag * rowColumnK->atU(p) ;
	 w->at(p) = diag * rowColumnK->atL(p) ;}

      // compute diagonal coefficient of rowColumn k
      rowColumnK->atDiag() -= rowColumnK->dot(r,'R',startK,k-1) ;
      diag = rowColumnK->atDiag() ;

      // test pivot not too small
      if (fabs(diag) < PRECISION) {
	 printf ("too small pivot %d\n",k) ;
	 exit(0) ;}

      // compute off-diagonal coefficients of rowColumns i>k
      for (i=k+1 ; i<=size ; i++) {
	 rowColumnI = this->giveRowColumn(i) ;
	 startI   = rowColumnI->giveStart() ;
	 if (startI <= k) {
	    start = max(startI,startK) ;
	    rowColumnI->atL(k) -= rowColumnI->dot(r,'R',start,k-1) ;
            rowColumnI->atL(k) /= diag ;
	    rowColumnI->atU(k) -= rowColumnI->dot(w,'C',start,k-1) ;
            rowColumnI->atU(k) /= diag ;}}

      delete r ;
      delete w ;}

   isFactorized = TRUE ;
   return this ;
}


Skyline*  SkylineUnsym :: forwardReductionWith (FloatArray* b)
   // Computes y, the solution of the system  L.y = b , where L is the
   // lower half of the receiver (assumed to be in a factorized form).
   // Returns the receiver. Note : y overwrites b.
{
   int        k,start ;
   RowColumn* rowColumnK ;

   if (! size)
      return this ;                                      // null size system
   if (b->giveSize() != size) 
      b -> growTo(size) ;

   for (k=1 ; k<=size ; k++) {
      rowColumnK  = this->giveRowColumn(k) ;
      start     = rowColumnK->giveStart() ;
      b->at(k) -= rowColumnK->dot(b,'R',start,k-1) ;}

   return this ;
}