349 if (!stateStorageInitialized_)
352 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
353 "LSQRIter::initialize(): Cannot initialize state storage!");
355 std::string errstr(
"LSQRIter::initialize(): Specified multivectors must have a consistent length and width.");
360 RCP<const MV> lhsMV = lp_->getLHS();
362 const bool debugSerialLSQR =
false;
364 if( debugSerialLSQR )
366 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> lhsNorm(1);
367 MVT::MvNorm( *lhsMV, lhsNorm );
368 std::cout <<
"initializeLSQR lhsNorm " << lhsNorm[0] << std::endl;
372 RCP<const MV> rhsMV = lp_->getInitPrecResVec();
373 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
374 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
375 MVT::MvAddMv( one, *rhsMV, zero, *U_, *U_);
377 RCP<const OP> M_left = lp_->getLeftPrec();
378 RCP<const OP> A = lp_->getOperator();
379 RCP<const OP> M_right = lp_->getRightPrec();
381 if( debugSerialLSQR )
383 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> rhsNorm(1);
384 MVT::MvNorm( *U_, rhsNorm );
385 std::cout <<
"initializeLSQR | U_ | : " << rhsNorm[0] << std::endl;
397 if ( M_left.is_null())
405 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
406 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
407 OPT::Apply (*A, *tempInRangeOfA, *V_,
CONJTRANS);
410 if (! M_right.is_null())
412 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
414 OPT::Apply (*M_right, *tempInDomainOfA, *V_,
CONJTRANS);
419 MVT::MvAddMv( one, *V_, zero, *V_, *W_);
421 frob_mat_norm_ = zero;
438 if (initialized_ ==
false) {
443 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
444 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
445 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
448 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> alpha(1);
449 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> beta(1);
450 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> xi(1);
453 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> wnorm2(1);
454 ScalarType rhobar, phibar, cs1, phi, rho, cs, sn, theta, xxnorm = MTzero, common;
455 ScalarType zetabar, sn1, psi, res = zero, ddnorm = zero, gamma, tau;
456 ScalarType anorm2 = zero;
457 ScalarType cs2 = -one, sn2 = zero, gammabar, zeta = zero, delta;
461 AV = MVT::Clone( *U_, 1);
462 Teuchos::RCP<MV> AtU;
463 AtU = MVT::Clone( *V_, 1);
464 const bool debugSerialLSQR =
false;
467 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
471 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1,
LSQRIterateFailure,
472 "LSQRIter::iterate(): current linear system has more than one vector!" );
476 MVT::MvNorm( *U_, beta );
477 if( SCT::real(beta[0]) > zero )
479 MVT::MvScale( *U_, one / beta[0] );
484 MVT::MvScale( *V_, one / beta[0] );
489 MVT::MvNorm( *V_, alpha );
490 if( debugSerialLSQR )
494 std::cout << iter_ <<
" First alpha " << alpha[0] <<
" beta " << beta[0] <<
" lambda " << lambda_ << std::endl;
496 if( SCT::real(alpha[0]) > zero )
498 MVT::MvScale( *V_, one / alpha[0] );
500 if( beta[0] * alpha[0] > zero )
502 MVT::MvScale( *W_, one / (beta[0] * alpha[0]) );
506 MVT::MvScale( *W_, zero );
510 RCP<const OP> M_left = lp_->getLeftPrec();
511 RCP<const OP> A = lp_->getOperator();
512 RCP<const OP> M_right = lp_->getRightPrec();
517 resid_norm_ = beta[0];
518 mat_resid_norm_ = alpha[0] * beta[0];
532 if ( M_right.is_null() )
534 OPT::Apply(*A, *V_, *AV);
538 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
539 OPT::Apply (*M_right, *V_, *tempInDomainOfA);
540 OPT::Apply(*A, *tempInDomainOfA, *AV);
543 if (! M_left.is_null())
545 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*AV);
546 OPT::Apply (*M_left, *tempInRangeOfA, *AV);
550 if ( !( M_left.is_null() && M_right.is_null() )
551 && debugSerialLSQR && iter_ == 1)
557 if (! M_left.is_null())
559 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
560 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
561 OPT::Apply (*A, *tempInRangeOfA, *AtU,
CONJTRANS);
567 if ( !( M_right.is_null() ) )
569 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
570 OPT::Apply (*M_right, *tempInDomainOfA, *AtU,
CONJTRANS);
573 MVT::MvAddMv( one, *AtU, -alpha[0], *V_, *AtU );
574 MVT::MvNorm( *AtU, xi );
575 std::cout <<
"| V alpha - A' u |= " << xi[0] << std::endl;
577 Teuchos::SerialDenseMatrix<int,ScalarType> uotuo(1,1);
578 MVT::MvTransMv( one, *U_, *U_, uotuo );
579 std::cout <<
"<U, U> = " << printMat(uotuo) << std::endl;
581 std::cout <<
"alpha = " << alpha[0] << std::endl;
583 Teuchos::SerialDenseMatrix<int,ScalarType> utav(1,1);
584 MVT::MvTransMv( one, *AV, *U_, utav );
585 std::cout <<
"<AV, U> = alpha = " << printMat(utav) << std::endl;
588 MVT::MvAddMv( one, *AV, -alpha[0], *U_, *U_ );
589 MVT::MvNorm( *U_, beta);
591 anorm2 += alpha[0]*alpha[0] + beta[0]*beta[0] + lambda_*lambda_;
594 if ( SCT::real(beta[0]) > zero )
597 MVT::MvScale( *U_, one / beta[0] );
599 if (M_left.is_null())
605 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
606 OPT::Apply (*M_left, *U_, *tempInRangeOfA,
CONJTRANS);
607 OPT::Apply(*A, *tempInRangeOfA, *AtU,
CONJTRANS);
609 if (! M_right.is_null())
611 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
612 OPT::Apply (*M_right, *tempInDomainOfA, *AtU,
CONJTRANS);
615 MVT::MvAddMv( one, *AtU, -beta[0], *V_, *V_ );
616 MVT::MvNorm( *V_, alpha );
623 if ( SCT::real(alpha[0]) > zero )
625 MVT::MvScale( *V_, one / alpha[0] );
630 common = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_);
631 cs1 = rhobar / common;
632 sn1 = lambda_ / common;
634 phibar = cs1 * phibar;
638 rho = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_ + beta[0]*beta[0]);
641 theta = sn * alpha[0];
642 rhobar = -cs * alpha[0];
644 phibar = sn * phibar;
648 gammabar = -cs2 * rho;
649 zetabar = (phi - delta*zeta) / gammabar;
650 sol_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(xxnorm + zetabar*zetabar);
651 gamma = Teuchos::ScalarTraits< ScalarType >::squareroot(gammabar*gammabar + theta*theta);
652 cs2 = gammabar / gamma;
654 zeta = (phi - delta*zeta) / gamma;
658 if ( M_right.is_null())
660 MVT::MvAddMv( phi / rho, *W_, one, *cur_soln_vec, *cur_soln_vec);
664 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*W_);
665 OPT::Apply (*M_right, *W_, *tempInDomainOfA);
666 MVT::MvAddMv( phi / rho, *tempInDomainOfA, one, *cur_soln_vec, *cur_soln_vec);
669 MVT::MvNorm( *W_, wnorm2 );
670 ddnorm += (one / rho)*(one / rho) * wnorm2[0]*wnorm2[0];
671 MVT::MvAddMv( one, *V_, -theta / rho, *W_, *W_ );
673 frob_mat_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(anorm2);
674 mat_cond_num_ = frob_mat_norm_ * Teuchos::ScalarTraits< ScalarType >::squareroot(ddnorm);
676 resid_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(phibar*phibar + res);
677 mat_resid_norm_ = alpha[0] * Teuchos::ScalarTraits< ScalarType >::magnitude(tau);