00001 #include "SundanceStochBlockJacobiSolver.hpp"
00002 #include "Sundance.hpp"
00003
00004
00005
00006 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00007 #include "PlayaVectorImpl.hpp"
00008 #endif
00009
00010
00011
00012 namespace Sundance
00013 {
00014
00015 void
00016 StochBlockJacobiSolver::solve(const Array<LinearOperator<double> >& KBlock,
00017 const Array<Vector<double> >& fBlock,
00018 Array<Vector<double> >& xBlock) const
00019 {
00020 Array<int> hasNonzeroMatrix(KBlock.size());
00021 for (int i=0; i<KBlock.size(); i++) hasNonzeroMatrix[i] = true;
00022
00023 solve(KBlock, hasNonzeroMatrix, fBlock, xBlock);
00024 }
00025
00026
00027 void
00028 StochBlockJacobiSolver::solve(const Array<LinearOperator<double> >& KBlock,
00029 const Array<int>& hasNonzeroMatrix,
00030 const Array<Vector<double> >& fBlock,
00031 Array<Vector<double> >& xBlock) const
00032 {
00033 int L = KBlock.size();
00034 int P = pcBasis_.nterms();
00035 int Q = fBlock.size();
00036
00037
00038
00039
00040 Array<Vector<double> > uPrev(P);
00041 Array<Vector<double> > uCur(P);
00042
00043 for (int i=0; i<P; i++)
00044 {
00045 TEUCHOS_TEST_FOR_EXCEPTION(fBlock[0].ptr().get()==0,
00046 std::runtime_error, "empty RHS vector block i=[" << i << "]");
00047 uPrev[i] = fBlock[0].copy();
00048 uCur[i] = fBlock[0].copy();
00049 uPrev[i].zero();
00050 uCur[i].zero();
00051 }
00052
00053 if (verbosity_) Out::root() << "starting Jacobi loop" << std::endl;
00054 bool converged = false;
00055 for (int iter=0; iter<maxIters_; iter++)
00056 {
00057 if (verbosity_) Out::root() << "Jacobi iter=" << iter << std::endl;
00058 bool haveNonConvergedBlock = false;
00059 double maxErr = -1.0;
00060 int numNonzeroBlocks = 0;
00061 for (int i=0; i<P; i++)
00062 {
00063 if (verbosity_) Out::root() << "Iter " << iter << ": block row i=" << i << " of " << P << " ..." << ends;
00064 Vector<double> b = fBlock[0].copy();
00065 b.zero();
00066 int nVecAdds = 0;
00067 for (int j=0; j<Q; j++)
00068 {
00069 double c_ij0 = pcBasis_.expectation(i,j,0);
00070 if (std::fabs(c_ij0) > 0.0)
00071 {
00072 b = b + c_ij0 * fBlock[j];
00073 nVecAdds++;
00074 }
00075 if (j>=L) continue;
00076 if (!hasNonzeroMatrix[j]) continue;
00077 Vector<double> tmp = fBlock[0].copy();
00078 tmp.zero();
00079 bool blockIsNeeded = false;
00080 for (int k=0; k<P; k++)
00081 {
00082 if (j==0 && k==i) continue;
00083 double c_ijk = pcBasis_.expectation(i,j,k);
00084 if (std::fabs(c_ijk) > 0.0)
00085 {
00086 tmp = tmp + c_ijk * uPrev[k];
00087 nVecAdds++;
00088 blockIsNeeded = true;
00089 }
00090 }
00091 numNonzeroBlocks += blockIsNeeded;
00092 b = (b - KBlock[j]*tmp);
00093 nVecAdds++;
00094 }
00095 b = b * (1.0/pcBasis_.expectation(i,i,0));
00096 if (verbosity_) Out::root() << "num vec adds = " << nVecAdds << std::endl;
00097 diagonalSolver_.solve(KBlock[0], b, uCur[i]);
00098 double err = (uCur[i]-uPrev[i]).norm2();
00099 if (err > convTol_) haveNonConvergedBlock=true;
00100 if (err > maxErr) maxErr = err;
00101 }
00102
00103
00104 for (int i=0; i<P; i++) uPrev[i] = uCur[i].copy();
00105
00106
00107 if (!haveNonConvergedBlock)
00108 {
00109 if (verbosity_) Out::root() << "=======> max err=" << maxErr << std::endl;
00110 if (verbosity_) Out::root() << "=======> converged! woo-hoo!" << std::endl;
00111 if (verbosity_) Out::root() << "estimated storage cost: "
00112 << setprecision(3) << 100*((double) L)/((double) numNonzeroBlocks)
00113 << " percent of monolithic storage" << std::endl;
00114 converged = true;
00115 break;
00116 }
00117 else
00118 {
00119 if (verbosity_) Out::root() << "maxErr=" << maxErr << ", trying again" << std::endl;
00120 }
00121 }
00122
00123 TEUCHOS_TEST_FOR_EXCEPT(!converged);
00124 xBlock = uCur;
00125 }
00126 }
