00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "LOCA_BorderedSolver_HouseholderQR.H"
00043 #include "LOCA_GlobalData.H"
00044 #include "LOCA_ErrorCheck.H"
00045
00046 LOCA::BorderedSolver::HouseholderQR::HouseholderQR(
00047 const Teuchos::RCP<LOCA::GlobalData>& global_data) :
00048 globalData(global_data),
00049 dblas()
00050 {
00051 }
00052
00053 LOCA::BorderedSolver::HouseholderQR::~HouseholderQR()
00054 {
00055 }
00056
00057 void
00058 LOCA::BorderedSolver::HouseholderQR::computeQR(
00059 const NOX::Abstract::MultiVector::DenseMatrix& C,
00060 const NOX::Abstract::MultiVector& B,
00061 bool use_c_transpose,
00062 NOX::Abstract::MultiVector::DenseMatrix& Y1,
00063 NOX::Abstract::MultiVector& Y2,
00064 NOX::Abstract::MultiVector::DenseMatrix& T,
00065 NOX::Abstract::MultiVector::DenseMatrix& R)
00066 {
00067 double beta;
00068 int m = B.numVectors();
00069
00070
00071 Y1.putScalar(0.0);
00072 T.putScalar(0.0);
00073 Y2 = B;
00074 if (use_c_transpose) {
00075 for (int i=0; i<m; i++)
00076 for (int j=0; j<m; j++)
00077 R(i,j) = C(j,i);
00078 }
00079 else
00080 R.assign(C);
00081
00082
00083 Teuchos::RCP<NOX::Abstract::MultiVector> v2 = Y2.clone(1);
00084
00085 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> v1;
00086 Teuchos::RCP<NOX::Abstract::MultiVector> h2;
00087 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> h1;
00088 Teuchos::RCP<NOX::Abstract::MultiVector> y2;
00089 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> y1;
00090 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> z;
00091 vector<int> h_idx;
00092 vector<int> y_idx;
00093 y_idx.reserve(m);
00094
00095 for (int i=0; i<m; i++) {
00096
00097
00098 v1 =
00099 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
00100 Y1,
00101 m-i,
00102 1, i, i));
00103
00104
00105 h_idx.resize(m-i);
00106 for (unsigned int j=0; j<h_idx.size(); j++)
00107 h_idx[j] = i+j;
00108 h2 = Y2.subView(h_idx);
00109
00110
00111 h1 =
00112 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
00113 R,
00114 m-i,
00115 m-i,
00116 i, i));
00117
00118 if (i > 0) {
00119
00120
00121 y_idx.push_back(i-1);
00122 y2 = Y2.subView(y_idx);
00123
00124
00125 y1 =
00126 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
00127 Y1,
00128 m-i,
00129 i, i, 0));
00130
00131
00132 z =
00133 Teuchos::rcp(new NOX::Abstract::MultiVector::DenseMatrix(Teuchos::View,
00134 T,
00135 i,
00136 1,
00137 0, i));
00138 }
00139
00140
00141 computeHouseholderVector(i, R, Y2, *v1, *v2, beta);
00142
00143
00144 applyHouseholderVector(*v1, *v2, beta, *h1, *h2);
00145
00146
00147 Y2[i] = (*v2)[0];
00148
00149 T(i,i) = -beta;
00150
00151 if (i > 0) {
00152
00153
00154 v2->multiply(1.0, *y2, *z);
00155
00156
00157 z->multiply(Teuchos::TRANS, Teuchos::NO_TRANS, -beta, *y1, *v1, -beta);
00158
00159
00160 dblas.TRMV(Teuchos::UPPER_TRI, Teuchos::NO_TRANS, Teuchos::NON_UNIT_DIAG,
00161 i, T.values(), m, z->values(), 1);
00162
00163 }
00164 }
00165
00166 }
00167
00168 void
00169 LOCA::BorderedSolver::HouseholderQR::computeHouseholderVector(
00170 int col,
00171 const NOX::Abstract::MultiVector::DenseMatrix& A1,
00172 const NOX::Abstract::MultiVector& A2,
00173 NOX::Abstract::MultiVector::DenseMatrix& V1,
00174 NOX::Abstract::MultiVector& V2,
00175 double& beta)
00176 {
00177 double houseP = A1(col,col);
00178
00179 V1(0,0) = 1.0;
00180 V2[0] = A2[col];
00181
00182 double sigma = A2[col].innerProduct(A2[col]);
00183 for (int i=col+1; i<A1.numRows(); i++)
00184 sigma += A1(i,col)*A1(i,col);
00185
00186 if (sigma == 0.0)
00187 beta = 0.0;
00188 else {
00189 double mu = sqrt(houseP*houseP + sigma);
00190 if (houseP <= 0.0)
00191 houseP = houseP - mu;
00192 else
00193 houseP = -sigma / (houseP + mu);
00194 beta = 2.0*houseP*houseP/(sigma + houseP*houseP);
00195
00196 V2.scale(1.0/houseP);
00197 for (int i=1; i<V1.numRows(); i++)
00198 V1(i,0) = A1(i+col,col) / houseP;
00199 }
00200
00201
00202 return;
00203 }
00204
00205 void
00206 LOCA::BorderedSolver::HouseholderQR::applyHouseholderVector(
00207 const NOX::Abstract::MultiVector::DenseMatrix& V1,
00208 const NOX::Abstract::MultiVector& V2,
00209 double beta,
00210 NOX::Abstract::MultiVector::DenseMatrix& A1,
00211 NOX::Abstract::MultiVector& A2)
00212 {
00213 int nColsA = A2.numVectors();
00214
00215
00216 NOX::Abstract::MultiVector::DenseMatrix u(1, nColsA);
00217 A2.multiply(1.0, V2, u);
00218
00219
00220 u.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, V1, A1, 1.0);
00221
00222
00223 A1.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -beta, V1, u, 1.0);
00224
00225
00226 A2.update(Teuchos::NO_TRANS, -beta, V2, u, 1.0);
00227 }
00228
00229 void
00230 LOCA::BorderedSolver::HouseholderQR::applyCompactWY(
00231 const NOX::Abstract::MultiVector::DenseMatrix& Y1,
00232 const NOX::Abstract::MultiVector& Y2,
00233 const NOX::Abstract::MultiVector::DenseMatrix& T,
00234 NOX::Abstract::MultiVector::DenseMatrix& X1,
00235 NOX::Abstract::MultiVector& X2,
00236 bool isZeroX1, bool isZeroX2,
00237 bool useTranspose) const
00238 {
00239 if (isZeroX1 && isZeroX2) {
00240 X1.putScalar(0.0);
00241 X2.init(0.0);
00242 return;
00243 }
00244
00245 int m = Y2.numVectors();
00246
00247 Teuchos::ETransp T_flag;
00248 if (useTranspose)
00249 T_flag = Teuchos::TRANS;
00250 else
00251 T_flag = Teuchos::NO_TRANS;
00252
00253 NOX::Abstract::MultiVector::DenseMatrix tmp(m, X2.numVectors());
00254
00255
00256 if (!isZeroX2)
00257 X2.multiply(1.0, Y2, tmp);
00258
00259
00260
00261 if (!isZeroX2 && !isZeroX1)
00262 tmp.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Y1, X1, 1.0);
00263 else if (!isZeroX1)
00264 tmp.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Y1, X1, 0.0);
00265
00266
00267 dblas.TRMM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, T_flag,
00268 Teuchos::NON_UNIT_DIAG, tmp.numRows(), tmp.numCols(), 1.0,
00269 T.values(), T.numRows(), tmp.values(), tmp.numRows());
00270
00271
00272
00273
00274 if (isZeroX1)
00275 X1.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Y1, tmp, 0.0);
00276 else
00277 X1.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Y1, tmp, 1.0);
00278
00279
00280 if (isZeroX2)
00281 X2.update(Teuchos::NO_TRANS, 1.0, Y2, tmp, 0.0);
00282 else
00283 X2.update(Teuchos::NO_TRANS, 1.0, Y2, tmp, 1.0);
00284 }
00285
00286 void
00287 LOCA::BorderedSolver::HouseholderQR::applyCompactWY(
00288 const NOX::Abstract::MultiVector::DenseMatrix& Y1,
00289 const NOX::Abstract::MultiVector& Y2,
00290 const NOX::Abstract::MultiVector::DenseMatrix& T,
00291 const NOX::Abstract::MultiVector::DenseMatrix* input1,
00292 const NOX::Abstract::MultiVector* input2,
00293 NOX::Abstract::MultiVector::DenseMatrix& result1,
00294 NOX::Abstract::MultiVector& result2,
00295 bool useTranspose) const
00296 {
00297 bool isZeroX1 = (input1 == NULL);
00298 bool isZeroX2 = (input2 == NULL);
00299
00300 if (!isZeroX1)
00301 result1.assign(*input1);
00302 if (!isZeroX2)
00303 result2 = *input2;
00304
00305 applyCompactWY(Y1, Y2, T, result1, result2, isZeroX1, isZeroX2,
00306 useTranspose);
00307 }
00308