83 finalP = MatrixFactory2::BuildCopy(rcpFromRef(P0));
88 "For now, solving Hessenberg system works only for a single iteration");
90 SC one = Teuchos::ScalarTraits<SC>::one(), zero = Teuchos::ScalarTraits<SC>::zero();
92 RCP<const Matrix> A = rcpFromRef(Aref);
99 ArrayRCP<const SC> D(A->getLocalNumRows(), one);
102 Teuchos::FancyOStream& mmfancy = this->GetOStream(
Statistics2);
105 RCP<Matrix> X = rcp_const_cast<Matrix>(rcpFromRef(P0)), tmpAP, newV;
106 std::vector<RCP<Matrix> > V(nIts_+1);
109 RCP<CrsMatrix> T_ = CrsMatrixFactory::Build(C.
GetPattern());
110 T_->fillComplete(P0.getDomainMap(), P0.getRangeMap());
111 RCP<Matrix> T = rcp(
new CrsMatrixWrap(T_));
117 tmpAP = MatrixMatrix::Multiply(*A,
false, *X,
false, mmfancy,
true,
true);
120 V[0] = MatrixFactory2::BuildCopy(T);
125 V[0]->scale(-one/rho);
128 std::vector<SC> h((nIts_+1) * (nIts_+1));
129 std::vector<SC> c(nIts_+1, 0.0);
130 std::vector<SC> s(nIts_+1, 0.0);
131 std::vector<SC> g(nIts_+1, 0.0);
134#define I(i,j) ((i) + (j)*(nIts_+1))
135 for (
size_t i = 0; i < nIts_; i++) {
137 tmpAP = MatrixMatrix::Multiply(*A,
false, *V[i],
false, mmfancy,
true,
true);
140 V[i+1] = MatrixFactory2::BuildCopy(T);
144 for (
size_t j = 0; j <= i; j++) {
148#ifndef TWO_ARG_MATRIX_ADD
149 newV = Teuchos::null;
150 MatrixMatrix::TwoMatrixAdd(*V[j],
false, -h[
I(j,i)], *V[i+1],
false, one, newV, mmfancy);
151 newV->fillComplete(V[i+1]->getDomainMap(), V[i+1]->getRangeMap());
162 MatrixMatrix::TwoMatrixAdd(*V[j],
false, -h[
I(j,i)], *V[i+1], one);
181 if (h[
I(i+1,i)] != zero) {
183 V[i+1]->scale(one/h[
I(i+1,i)]);
187 givapp(&c[0], &s[0], &h[
I(0,i)], i);
189 SC nu = sqrt(h[
I(i,i)]*h[
I(i,i)] + h[
I(i+1,i)]*h[
I(i+1,i)]);
191 c[i] = h[
I(i, i)] / nu;
192 s[i] = -h[
I(i+1,i)] / nu;
193 h[
I(i,i)] = c[i] * h[
I(i,i)] - s[i] * h[
I(i+1,i)];
196 givapp(&c[i], &s[i], &g[i], 1);
202 std::vector<SC> y(nIts_);
204 y[0] = g[0] / h[
I(0,0)];
209 for (
size_t i = 0; i < nIts_; i++) {
210#ifndef TWO_ARG_MATRIX_ADD
211 newV = Teuchos::null;
212 MatrixMatrix::TwoMatrixAdd(*V[i],
false, y[i], *finalP,
false, one, newV, mmfancy);
213 newV->fillComplete(finalP->getDomainMap(), finalP->getRangeMap());
216 MatrixMatrix::TwoMatrixAdd(*V[i],
false, y[i], *finalP, one);
static void MyOldScaleMatrix(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< const Scalar > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)