193 const value_type& alpha,
196 const value_type& beta)
198#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
199 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::GMRESDivisionStrategy::divide()");
202 ordinal_type sz = basis->size();
203 ordinal_type pa = a.
size();
204 ordinal_type pb = b.
size();
213 const value_type* ca = a.
coeff();
214 const value_type* cb = b.
coeff();
215 value_type* cc = c.
coeff();
222 if (pb < Cijk->num_k())
223 k_end = Cijk->find_k(pb);
229 j_it != Cijk->j_end(k_it); ++j_it) {
232 i_it != Cijk->i_end(j_it); ++i_it) {
235 (*A)(i,
j) += cijk*cb[k];
242 for (ordinal_type i=0; i<pa; i++)
243 (*B)(i,0) = ca[i]*basis->norm_squared(i);
245 Teuchos::SerialDenseMatrix<ordinal_type,value_type> D(sz, 1);
249 for (ordinal_type i=0; i<sz; i++){
250 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r(Teuchos::View, *A, 1, sz, i, 0);
251 D(i,0)=
sqrt(r.normOne());
254 for (ordinal_type i=0; i<sz; i++){
255 for (ordinal_type
j=0;
j<sz;
j++){
256 (*A)(i,
j)=(*A)(i,
j)/(D(i,0)*D(
j,0));
261 for (ordinal_type i=0; i<sz; i++){
262 (*B)(i,0)=(*B)(i,0)/D(i,0);
269 pb = basis->dimension()+1;
271 if (pb < Cijk->num_k())
272 k_end = Cijk->find_k(pb);
276 j_it != Cijk->j_end(k_it); ++j_it) {
279 i_it != Cijk->i_end(j_it); ++i_it) {
282 (*M)(i,
j) += cijk*cb[k];
290 for (ordinal_type i=0; i<sz; i++){
291 for (ordinal_type
j=0;
j<sz;
j++){
292 (*M)(i,
j)=(*M)(i,
j)/(D(i,0)*D(
j,0));
299 GMRES(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *M,
diag);
303 GMRES(*A,*X,*B, max_it, tol, prec_iter, basis->order(), basis->dimension(), PrecNum, *A,
diag);
308 for (ordinal_type i=0; i<sz; i++){
309 (*X)(i,0)=(*X)(i,0)/D(i,0);
314 for (ordinal_type i=0; i<pc; i++)
315 cc[i] = alpha*(*X)(i,0) + beta*cc[i];
318 for (ordinal_type i=0; i<pc; i++)
319 cc[i] = alpha*ca[i]/cb[0] + beta*cc[i];
327GMRES(
const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & A,
328 Teuchos::SerialDenseMatrix<ordinal_type,value_type> & X,
329 const Teuchos::SerialDenseMatrix<ordinal_type,value_type> & B,
330 ordinal_type max_iter,
331 value_type tolerance,
332 ordinal_type prec_iter,
335 ordinal_type PrecNum,
336 const Teuchos::SerialDenseMatrix<ordinal_type, value_type> & M,
339 ordinal_type n = A.numRows();
342 Teuchos::SerialDenseMatrix<ordinal_type, value_type> P(n,n);
343 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Ax(n,1);
344 Ax.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,1.0, A, X, 0.0);
345 Teuchos::SerialDenseMatrix<ordinal_type, value_type> r0(Teuchos::Copy,B);
347 resid=r0.normFrobenius();
349 Teuchos::SerialDenseMatrix<ordinal_type, value_type> v(n,1);
351 Teuchos::SerialDenseMatrix<ordinal_type, value_type> h(1,1);
353 Teuchos::SerialDenseMatrix<ordinal_type, value_type> V(n,1);
355 for (ordinal_type i=0; i<n; i++){
359 Teuchos::SerialDenseMatrix<ordinal_type, value_type> bb(1,1);
361 Teuchos::SerialDenseMatrix<ordinal_type, value_type> w(n,1);
362 Teuchos::SerialDenseMatrix<ordinal_type, value_type> c;
363 Teuchos::SerialDenseMatrix<ordinal_type, value_type> s;
364 while (resid > tolerance && k < max_iter){
369 Teuchos::SerialDenseMatrix<ordinal_type, value_type> vk(Teuchos::Copy, V, n,1,0,k-1);
371 Teuchos::SerialDenseMatrix<ordinal_type, value_type> z(vk);
376 else if (PrecNum == 2){
380 else if (PrecNum == 3){
384 else if (PrecNum == 4){
389 w.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1, A, z, 0.0);
390 Teuchos::SerialDenseMatrix<ordinal_type, value_type> vi(n,1);
391 Teuchos::SerialDenseMatrix<ordinal_type, value_type> ip(1,1);
392 for (ordinal_type i=0; i<k; i++){
394 Teuchos::SerialDenseMatrix<ordinal_type, value_type> vi(Teuchos::Copy, V, n,1,0,i);
396 ip.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, vi, w, 0.0);
402 h(k,k-1)=w.normFrobenius();
403 w.scale(1.0/h(k,k-1));
405 for (ordinal_type i=0; i<n; i++){
410 for (ordinal_type i=0; i<k-1; i++){
411 value_type q=c(i,0)*h(i,k-1)+s(i,0)*h(i+1,k-1);
412 h(i+1,k-1)=-1*s(i,0)*h(i,k-1)+c(i,0)*h(i+1,k-1);
420 value_type l =
sqrt(h(k-1,k-1)*h(k-1,k-1)+h(k,k-1)*h(k,k-1));
421 c(k-1,0)=h(k-1,k-1)/l;
428 bb(k,0)=-s(k-1,0)*bb(k-1,0);
429 bb(k-1,0)=c(k-1,0)*bb(k-1,0);
432 resid =
fabs(bb(k,0));
436 bb.reshape(h.numRows()-1 ,1);
440 lapack.TRTRS(
'U',
'N',
'N', h.numRows()-1, 1, h.values(), h.stride(), bb.values(), bb.stride(),&info);
441 Teuchos::SerialDenseMatrix<ordinal_type, value_type> ans(X);
443 ans.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, V, bb, 0.0);
448 else if (PrecNum == 2){
452 else if (PrecNum == 3){
456 else if (PrecNum == 4){
ordinal_type GMRES(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &X, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &B, ordinal_type max_iter, value_type tolerance, ordinal_type prec_iter, ordinal_type order, ordinal_type dim, ordinal_type PrecNum, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &M, ordinal_type diag)