312 Time_->ResetStartTime();
318 if (PolyDegree_ <= 0)
321#ifdef HAVE_IFPACK_EPETRAEXT
323 if(IsRowMatrix_ && InvDiagonal_ == Teuchos::null && UseBlockMode_){
327 if(!CrsMatrix) UseBlockMode_=
false;
330 InvBlockDiagonal_=Teuchos::rcp(
new EpetraExt_PointToBlockDiagPermute(*CrsMatrix));
331 if(InvBlockDiagonal_==Teuchos::null) IFPACK_CHK_ERR(-6);
333 ierr=InvBlockDiagonal_->SetParameters(BlockList_);
334 if(ierr) IFPACK_CHK_ERR(-7);
336 ierr=InvBlockDiagonal_->Compute();
337 if(ierr) IFPACK_CHK_ERR(-8);
341 if (IsRowMatrix_ && InvDiagonal_ == Teuchos::null && !UseBlockMode_)
345 if (InvDiagonal_ == Teuchos::null)
348 IFPACK_CHK_ERR(
Matrix().ExtractDiagonalCopy(*InvDiagonal_));
352 for (
int i = 0 ; i < NumMyRows_ ; ++i) {
353 double diag = (*InvDiagonal_)[i];
354 if (IFPACK_ABS(diag) < MinDiagonalValue_)
355 (*InvDiagonal_)[i] = MinDiagonalValue_;
357 (*InvDiagonal_)[i] = 1.0 / diag;
362 double lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max;
363 if (LambdaRealMax_ == -1) {
365 GMRES(
Matrix(), *InvDiagonal_, EigMaxIters_, lambda_real_min, lambda_real_max, lambda_imag_min, lambda_imag_max);
366 LambdaRealMin_=lambda_real_min; LambdaImagMin_=lambda_imag_min;
367 LambdaRealMax_=lambda_real_max; LambdaImagMax_=lambda_imag_max;
378 const std::complex<double> zero(0.0,0.0);
381 double lenx = LambdaRealMax_-LambdaRealMin_;
382 int nx = ceil(lenx*((
double) LSPointsReal_));
383 if (nx<2) { nx = 2; }
384 double hx = lenx/((double) nx);
385 std::vector<double> xs;
386 if(abs(lenx)>1.0e-8) {
387 for(
int pt=0; pt<=nx; pt++ ) {
388 xs.push_back(hx*pt+LambdaRealMin_);
392 xs.push_back(LambdaRealMax_);
395 double leny = LambdaImagMax_-LambdaImagMin_;
396 int ny = ceil(leny*((
double) LSPointsImag_));
397 if (ny<2) { ny = 2; }
398 double hy = leny/((double) ny);
399 std::vector<double> ys;
400 if(abs(leny)>1.0e-8) {
401 for(
int pt=0; pt<=ny; pt++ ) {
402 ys.push_back(hy*pt+LambdaImagMin_);
406 ys.push_back(LambdaImagMax_);
409 std::vector< std::complex<double> > cpts;
410 for(
int jj=0; jj<ny; jj++ ) {
411 for(
int ii=0; ii<nx; ii++ ) {
412 std::complex<double> cpt(xs[ii],ys[jj]);
416 cpts.push_back(zero);
418#ifdef HAVE_TEUCHOS_COMPLEX
419 const std::complex<double> one(1.0,0.0);
422 Teuchos::SerialDenseMatrix<int, std::complex<double> > Vmatrix(cpts.size(),PolyDegree_+1);
423 Vmatrix.putScalar(zero);
424 for (
int jj = 0; jj <= PolyDegree_; ++jj) {
425 for (
int ii = 0; ii < static_cast<int> (cpts.size ()) - 1; ++ii) {
427 Vmatrix(ii,jj) = pow(cpts[ii],jj);
430 Vmatrix(ii,jj) = one;
434 Vmatrix(cpts.size()-1,0)=one;
437 Teuchos::SerialDenseMatrix< int,std::complex<double> > RHS(cpts.size(),1);
439 RHS(cpts.size()-1,0)=one;
442 Teuchos::LAPACK< int, std::complex<double> > lapack;
443 const int N = Vmatrix.numCols();
444 Teuchos::Array<double> singularValues(N);
445 Teuchos::Array<double> rwork(1);
446 rwork.resize (std::max (1, 5 * N));
447 std::complex<double> lworkScalar(1.0,0.0);
449 lapack.GELS(
'N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
450 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
451 &lworkScalar, -1, &info);
452 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
453 "_GELSS workspace query returned INFO = "
454 << info <<
" != 0.");
455 const int lwork =
static_cast<int> (real(lworkScalar));
456 TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
457 "_GELSS workspace query returned LWORK = "
458 << lwork <<
" < 0.");
460 Teuchos::Array< std::complex<double> > work (std::max (1, lwork));
462 lapack.GELS(
'N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
463 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
464 &work[0], lwork, &info);
466 coeff_.resize(PolyDegree_+1);
467 std::complex<double> c0=RHS(0,0);
468 for(
int ii=0; ii<=PolyDegree_; ii++) {
473 coeff_[ii]=real(RHS(ii,0)/c0);
480 Teuchos::SerialDenseMatrix< int, double > Vmatrix(xs.size()+1,PolyDegree_+1);
481 Vmatrix.putScalar(0.0);
482 for(
int jj=0; jj<=PolyDegree_; jj++) {
483 for( std::vector<double>::size_type ii=0; ii<xs.size(); ii++) {
485 Vmatrix(ii,jj)=pow(xs[ii],jj);
492 Vmatrix(xs.size(),0)=1.0;
495 Teuchos::SerialDenseMatrix< int, double > RHS(xs.size()+1,1);
497 RHS(xs.size(),0)=1.0;
500 Teuchos::LAPACK< int, double > lapack;
501 const int N = Vmatrix.numCols();
502 Teuchos::Array<double> singularValues(N);
503 Teuchos::Array<double> rwork(1);
504 rwork.resize (std::max (1, 5 * N));
505 double lworkScalar(1.0);
507 lapack.GELS(
'N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
508 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
509 &lworkScalar, -1, &info);
510 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
511 "_GELSS workspace query returned INFO = "
512 << info <<
" != 0.");
513 const int lwork =
static_cast<int> (lworkScalar);
514 TEUCHOS_TEST_FOR_EXCEPTION(lwork < 0, std::logic_error,
515 "_GELSS workspace query returned LWORK = "
516 << lwork <<
" < 0.");
518 Teuchos::Array< double > work (std::max (1, lwork));
520 lapack.GELS(
'N', Vmatrix.numRows(), Vmatrix.numCols(), RHS.numCols(),
521 Vmatrix.values(), Vmatrix.numRows(), RHS.values(), RHS.numRows(),
522 &work[0], lwork, &info);
524 coeff_.resize(PolyDegree_+1);
526 for(
int ii=0; ii<=PolyDegree_; ii++) {
531 coeff_[ii]=RHS(ii,0)/c0;
536#ifdef IFPACK_FLOPCOUNTERS
537 ComputeFlops_ += NumMyRows_;
541 ComputeTime_ += Time_->ElapsedTime();