← Back to branch summary

~chaffra/+junk/trilinos

« back to all changes in this revision

Viewing changes to packages/anasazi/src/AnasaziSIRTR.hpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Johannes Ring
  • Date: 2009-12-13 12:53:22 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20091213125322-in0nrdjc55deqsw9
Tags: 10.0.3.dfsg-1
[Christophe Prud'homme]
* New upstream release

[Johannes Ring]
* debian/patches/libname.patch: Add prefix 'libtrilinos_' to all
  libraries. 
* debian/patches/soname.patch: Add soversion to libraries.
* debian/watch: Update download URL.
* debian/control:
  - Remove python-numeric from Build-Depends (virtual package).
  - Remove automake and autotools from Build-Depends and add cmake to
    reflect switch to CMake.
  - Add python-support to Build-Depends.
* debian/rules: 
  - Cleanup and updates for switch to CMake.

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/anasazi/src/AnasaziSIRTR.hpp
130
130
    trRetType innerStop_;
131
131
    // 
132
132
    // number of inner iterations
133
 
    int innerIters_;
 
133
    int innerIters_, totalInnerIters_;
134
134
  };
135
135
 
136
136
 
149
149
        ) : 
150
150
    RTRBase<ScalarType,MV,OP>(problem,sorter,printer,tester,ortho,params,"SIRTR",true), 
151
151
    ZERO(MAT::zero()),
152
 
    ONE(MAT::one())
 
152
    ONE(MAT::one()),
 
153
    totalInnerIters_(0)
153
154
  {     
154
155
    // set up array of stop reasons
155
156
    stopReasons_.push_back("n/a");
262
263
    }
263
264
 
264
265
    // 
265
 
    // The preconditioner is 
 
266
    // For Olsen preconditioning, the preconditioner is 
266
267
    // Z = P_{Prec^-1 BX, BX} Prec^-1 R
267
268
    // for efficiency, we compute Prec^-1 BX once here for use later
268
 
    if (this->hasPrec_) 
 
269
    // Otherwise, we don't need PBX
 
270
    if (this->hasPrec_ && this->olsenPrec_) 
269
271
    {
270
272
      std::vector<int> ind(this->blockSize_);
271
273
      for (int i=0; i<this->blockSize_; ++i) ind[i] = this->numAuxVecs_+i;
279
281
    }
280
282
 
281
283
    // Z = P_{Prec^-1 BV, BV} Prec^-1 R
282
 
    // Prec^-1 BV in PBV
 
284
    //    Prec^-1 BV in PBV
 
285
    // or
 
286
    // Z = P_{BV,BV} Prec^-1 R
283
287
    if (this->hasPrec_) 
284
288
    {
285
289
      TimeMonitor prectimer( *this->timerPrec_ );
286
290
      OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
287
291
      this->counterPrec_ += this->blockSize_;
288
292
      // the orthogonalization time counts under Ortho and under Preconditioning
289
 
      {
 
293
      if (this->olsenPrec_) {
290
294
        TimeMonitor orthtimer( *this->timerOrtho_ );
291
295
        this->orthman_->projectGen(
292
296
            *this->Z_,                                             // operating on Z
294
298
            tuple<PSDM>(null),                                     // don't care about coeffs
295
299
            null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));       // don't have B*PBV or B*Z, but do have B*V
296
300
      }
 
301
      else {
 
302
        TimeMonitor orthtimer( *this->timerOrtho_ );
 
303
        this->orthman_->projectGen(
 
304
            *this->Z_,                                             // operating on Z
 
305
            tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,    // P_{BV,V}, and <BV,V>_B != I
 
306
            tuple<PSDM>(null),                                     // don't care about coeffs
 
307
            null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));       // don't have B*BV, but do have B*V
 
308
      }
297
309
      ginnersep(*this->R_,*this->Z_,z_r);
298
310
    }
299
311
    else {
607
619
      }
608
620
 
609
621
      // Z = P_{Prec^-1 BV, BV} Prec^-1 R
610
 
      // Prec^-1 BV in PBV
 
622
      //    Prec^-1 BV in PBV
 
623
      // or
 
624
      // Z = P_{BV,BV} Prec^-1 R
611
625
      zold_rold = z_r;
612
626
      if (this->hasPrec_)
613
627
      {
615
629
        OPT::Apply(*this->Prec_,*this->R_,*this->Z_);
616
630
        this->counterPrec_ += this->blockSize_;
617
631
        // the orthogonalization time counts under Ortho and under Preconditioning
618
 
        {
 
632
        if (this->olsenPrec_) {
619
633
          TimeMonitor orthtimer( *this->timerOrtho_ );
620
634
          this->orthman_->projectGen(
621
635
              *this->Z_,                                             // operating on Z
623
637
              tuple<PSDM>(null),                                     // don't care about coeffs
624
638
              null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));       // don't have B*PBV or B*Z, but do have B*V
625
639
        }
 
640
        else {
 
641
          TimeMonitor orthtimer( *this->timerOrtho_ );
 
642
          this->orthman_->projectGen(
 
643
              *this->Z_,                                             // operating on Z
 
644
              tuple<PCMV>(this->BV_),tuple<PCMV>(this->V_),false,    // P_{BV,V}, and <BV,V>_B != I
 
645
              tuple<PSDM>(null),                                     // don't care about coeffs
 
646
              null,tuple<PCMV>(null), tuple<PCMV>(this->BV_));       // don't have B*BV, but do have B*V
 
647
        }
626
648
        ginnersep(*this->R_,*this->Z_,z_r);
627
649
      }
628
650
      else {
734
756
 
735
757
      // solve the trust-region subproblem
736
758
      solveTRSubproblem();
 
759
      totalInnerIters_ += innerIters_;
737
760
 
738
761
      // perform debugging on eta et al.
739
762
      if (this->om_->isVerbosity( Debug ) ) {
762
785
 
763
786
 
764
787
      // compute the retraction of eta: R_X(eta) = X+eta
765
 
      // we must accept, but we will work out of delta so that we can multiply back into X below
 
788
      // we must accept, but we will work out of temporary so that we can multiply back into X below
766
789
      RCP<MV> XpEta;
767
790
      SIRTR_GET_TEMP_MV(XpEta,workspace);                 // workspace size is 3
768
791
      {
771
794
      }
772
795
 
773
796
      //
774
 
      // perform rayleigh-ritz for newX = X+eta
 
797
      // perform rayleigh-ritz for XpEta = X+eta
775
798
      // save an old copy of f(X) for rho analysis below
776
799
      //
777
800
      MagnitudeType oldfx = this->fx_;
808
831
        }
809
832
      }
810
833
      else {
811
 
        // project B onto X+eta
 
834
        // project I onto X+eta
812
835
        TimeMonitor lcltimer( *this->timerLocalProj_ );
813
836
        MVT::MvTransMv(ONE,*XpEta,*XpEta,BB);
814
837
      }
822
845
        ret = Utils::directSolver(this->blockSize_,AA,Teuchos::rcpFromRef(BB),S,this->theta_,rank,1);
823
846
      }
824
847
      this->om_->stream(Debug) << "S: " << std::endl << S << std::endl;;
825
 
      TEST_FOR_EXCEPTION(ret != 0,std::logic_error,"Anasazi::SIRTR::iterate(): failure solving projected eigenproblem after retraction. ret == " << ret);
 
848
      TEST_FOR_EXCEPTION(ret != 0,std::logic_error,"Anasazi::SIRTR::iterate(): failure solving projected eigenproblem after retraction. ret == " << ret << "AA: " << AA << std::endl << "BB: " << BB << std::endl);
826
849
      TEST_FOR_EXCEPTION(rank != this->blockSize_,RTRRitzFailure,"Anasazi::SIRTR::iterate(): retracted iterate failed in Ritz analysis. rank == " << rank);
827
850
 
828
851
      //
1028
1051
    os <<"Parameter rho_prime is  " << rho_prime_ << endl;
1029
1052
    os <<"Inner stopping condition was " << stopReasons_[innerStop_] << endl;
1030
1053
    os <<"Number of inner iterations was " << innerIters_ << endl;
 
1054
    os <<"Total number of inner iterations is " << totalInnerIters_ << endl;
1031
1055
    os <<"f(x) is " << this->fx_ << endl;
1032
1056
 
1033
1057
    os.setf(std::ios_base::right, std::ios_base::adjustfield);