C++ LAPACK::LAPY2方法代码示例

NOX::Abstract::Group::ReturnType
LOCA::EigenvalueSort::LargestMagnitude::sort(int n, double* r_evals,
                         double* i_evals,
                         std::vector<int>* perm) const
{
  int i, j;
  int tempord = 0;
  double temp, tempr, tempi;
  Teuchos::LAPACK<int,double> lapack;

  //
  // Reset the index
  if (perm) {
    for (i=0; i < n; i++) {
      (*perm)[i] = i;
    }
  }

  //---------------------------------------------------------------
  // Sort eigenvalues in decreasing order of magnitude
  //---------------------------------------------------------------
  for (j=1; j < n; ++j) {
    tempr = r_evals[j]; tempi = i_evals[j];
    if (perm)
      tempord = (*perm)[j];
    temp=lapack.LAPY2(r_evals[j],i_evals[j]);
    for (i=j-1; i>=0 && lapack.LAPY2(r_evals[i],i_evals[i])<temp; --i) {
      r_evals[i+1]=r_evals[i]; i_evals[i+1]=i_evals[i];
      if (perm)
    (*perm)[i+1]=(*perm)[i];
    }
    r_evals[i+1] = tempr; i_evals[i+1] = tempi;
    if (perm)
      (*perm)[i+1] = tempord;
  }
  return NOX::Abstract::Group::Ok;
}

//.........这里部分代码省略.........
      std::vector<int> curind(1);
      std::vector<double> resnorm(1), tempnrm(1);
      Teuchos::RCP<MV> tempAevec;
      Teuchos::RCP<const MV> evecr, eveci;
      Epetra_MultiVector Aevec(Map,numev);
      
      // Compute A*evecs
      OPT::Apply( *A, *evecs, Aevec );
      
      Teuchos::SerialDenseMatrix<int,double> Breal(1,1), Bimag(1,1);
      while (i<numev) {
        if (index[i]==0) {
          // Get a view of the current eigenvector (evecr)
          curind[0] = i;
          evecr = MVT::CloneView( *evecs, curind );
	  
          // Get a copy of A*evecr
          tempAevec = MVT::CloneCopy( Aevec, curind );
	  
          // Compute A*evecr - lambda*evecr
          Breal(0,0) = evals[i].realpart;
          MVT::MvTimesMatAddMv( -1.0, *evecr, Breal, 1.0, *tempAevec );
	  
          // Compute the norm of the residual and increment counter
          MVT::MvNorm( *tempAevec, resnorm );
          normA[i] = resnorm[0]/Teuchos::ScalarTraits<MagnitudeType>::magnitude( evals[i].realpart );
          i++;
        } else {
          // Get a view of the real part of the eigenvector (evecr)
          curind[0] = i;
          evecr = MVT::CloneView( *evecs, curind );
	  
          // Get a copy of A*evecr
          tempAevec = MVT::CloneCopy( Aevec, curind );
	  
          // Get a view of the imaginary part of the eigenvector (eveci)
          curind[0] = i+1;
          eveci = MVT::CloneView( *evecs, curind );
	  
          // Set the eigenvalue into Breal and Bimag
          Breal(0,0) = evals[i].realpart;
          Bimag(0,0) = evals[i].imagpart;
	  
          // Compute A*evecr - evecr*lambdar + eveci*lambdai
          MVT::MvTimesMatAddMv( -1.0, *evecr, Breal, 1.0, *tempAevec );
          MVT::MvTimesMatAddMv( 1.0, *eveci, Bimag, 1.0, *tempAevec );
          MVT::MvNorm( *tempAevec, tempnrm );
	  
          // Get a copy of A*eveci
          tempAevec = MVT::CloneCopy( Aevec, curind );
	  
          // Compute A*eveci - eveci*lambdar - evecr*lambdai
          MVT::MvTimesMatAddMv( -1.0, *evecr, Bimag, 1.0, *tempAevec );
          MVT::MvTimesMatAddMv( -1.0, *eveci, Breal, 1.0, *tempAevec );
          MVT::MvNorm( *tempAevec, resnorm );
	  
          // Compute the norms and scale by magnitude of eigenvalue
          normA[i] = lapack.LAPY2( tempnrm[i], resnorm[i] ) /
            lapack.LAPY2( evals[i].realpart, evals[i].imagpart );
          normA[i+1] = normA[i];
	  
          i=i+2;
        }
      }
    }
    
    // Output computed eigenvalues and their direct residuals
    if (verbose && MyPID==0) {
      cout.setf(std::ios_base::right, std::ios_base::adjustfield);
      cout<<endl<< "Actual Residuals"<<endl;
      if (MyProblem->isHermitian()) {
        cout<< std::setw(16) << "Real Part"
	    << std::setw(20) << "Direct Residual"<< endl;
        cout<<"-----------------------------------------------------------"<<endl;
        for (int i=0; i<numev; i++) {
          cout<< std::setw(16) << evals[i].realpart 
	      << std::setw(20) << normA[i] << endl;
        }  
        cout<<"-----------------------------------------------------------"<<endl;
      } 
      else {
        cout<< std::setw(16) << "Real Part"
	    << std::setw(16) << "Imag Part"
	    << std::setw(20) << "Direct Residual"<< endl;
        cout<<"-----------------------------------------------------------"<<endl;
        for (int i=0; i<numev; i++) {
          cout<< std::setw(16) << evals[i].realpart 
	      << std::setw(16) << evals[i].imagpart 
	      << std::setw(20) << normA[i] << endl;
        }  
        cout<<"-----------------------------------------------------------"<<endl;
      }  
    }
  }
  
#ifdef EPETRA_MPI
  MPI_Finalize();
#endif  
  return 0;
}

//.........这里部分代码省略.........
    Epetra_MultiVector Aevec(Map,numev);

    // Compute A*evecs
    OpTraits::Apply( *A, *evecs, Aevec );

    Teuchos::SerialDenseMatrix<int,double> Breal(1,1), Bimag(1,1);
    while (i<numev) {
      if (index[i]==0) {
        // Get a view of the current eigenvector (evecr)
        curind[0] = i;
        evecr = MVTraits::CloneView( *evecs, curind );

        // Get a copy of A*evecr
        tempAevec = MVTraits::CloneCopy( Aevec, curind );

        // Compute A*evecr - lambda*evecr
        Breal(0,0) = evals[i].realpart;
        MVTraits::MvTimesMatAddMv( -1.0, *evecr, Breal, 1.0, *tempAevec );

        // Compute the norm of the residual and increment counter
        MVTraits::MvNorm( *tempAevec, resnorm );
        normA[i] = resnorm[0] / Teuchos::ScalarTraits<MagnitudeType>::magnitude( evals[i].realpart );
        i++;
      } else {
        // Get a view of the real part of the eigenvector (evecr)
        curind[0] = i;
        evecr = MVTraits::CloneView( *evecs, curind );

        // Get a copy of A*evecr
        tempAevec = MVTraits::CloneCopy( Aevec, curind );

        // Get a view of the imaginary part of the eigenvector (eveci)
        curind[0] = i+1;
        eveci = MVTraits::CloneView( *evecs, curind );

        // Set the eigenvalue into Breal and Bimag
        Breal(0,0) = evals[i].realpart;
        Bimag(0,0) = evals[i].imagpart;

        // Compute A*evecr - evecr*lambdar + eveci*lambdai
        MVTraits::MvTimesMatAddMv( -1.0, *evecr, Breal, 1.0, *tempAevec );
        MVTraits::MvTimesMatAddMv( 1.0, *eveci, Bimag, 1.0, *tempAevec );
        MVTraits::MvNorm( *tempAevec, tempnrm );

        // Get a copy of A*eveci
        tempAevec = MVTraits::CloneCopy( Aevec, curind );

        // Compute A*eveci - eveci*lambdar - evecr*lambdai
        MVTraits::MvTimesMatAddMv( -1.0, *evecr, Bimag, 1.0, *tempAevec );
        MVTraits::MvTimesMatAddMv( -1.0, *eveci, Breal, 1.0, *tempAevec );
        MVTraits::MvNorm( *tempAevec, resnorm );

        // Compute the norms and scale by magnitude of eigenvalue
        normA[i] = lapack.LAPY2( tempnrm[0], resnorm[0] ) /
          lapack.LAPY2( evals[i].realpart, evals[i].imagpart );
        normA[i+1] = normA[i];

        i=i+2;
      }
    }

    // Output computed eigenvalues and their direct residuals
    if (verbose && MyPID==0) {
      std::cout.setf(std::ios_base::right, std::ios_base::adjustfield);
      std::cout<<std::endl<< "Actual Residuals"<<std::endl;
      std::cout<< std::setw(16) << "Real Part"
          << std::setw(16) << "Imag Part"
          << std::setw(20) << "Direct Residual"<< std::endl;
      std::cout<<"-----------------------------------------------------------"<<std::endl;
      for (int j=0; j<numev; j++) {
        std::cout<< std::setw(16) << evals[j].realpart
            << std::setw(16) << evals[j].imagpart
            << std::setw(20) << normA[j] << std::endl;
        if ( normA[j] > tol ) {
          testFailed = true;
        }
      }
      std::cout<<"-----------------------------------------------------------"<<std::endl;
    }
  }

#ifdef EPETRA_MPI
  MPI_Finalize();
#endif

  if (testFailed) {
    if (verbose && MyPID==0) {
      std::cout << "End Result: TEST FAILED" << std::endl;
    }
    return -1;
  }
  //
  // Default return value
  //
  if (verbose && MyPID==0) {
    std::cout << "End Result: TEST PASSED" << std::endl;
  }

  return 0;
}

//.........这里部分代码省略.........
      MyProblem->setHermitian(false);
      MyProblem->setNEV( nev );
      // Inform the eigenproblem that you are finished passing it information
      TEUCHOS_TEST_FOR_EXCEPTION( MyProblem->setProblem() != true, 
                          std::runtime_error, "Anasazi::BasicEigenproblem::setProblem() returned with error.");

      // Initialize the Block Arnoldi solver
      Anasazi::BlockKrylovSchurSolMgr<ST,MV,OP> MySolverMgr(MyProblem, MyPL);

      // Solve the problem to the specified tolerances or length
      Anasazi::ReturnType returnCode = MySolverMgr.solve();

      // Get the eigenvalues and eigenvectors from the eigenproblem
      Anasazi::Eigensolution<ST,MV> sol = MyProblem->getSolution();
      vector<Anasazi::Value<ST> > evals = sol.Evals;
      RCP<MV> evecs = sol.Evecs;
      vector<int> index = sol.index;
      int numev = sol.numVecs;

      if (returnCode != Anasazi::Converged) {
        cout << "Anasazi::EigensolverMgr::solve() returned unconverged, computing " << numev << " of " << nev << endl;
      }
      else {
        cout << "Anasazi::EigensolverMgr::solve() returned converged, computing " << numev << " of " << nev << endl;
      }

      if (numev > 0) {
        // Compute residuals.
        Teuchos::LAPACK<int,MT> lapack;
        vector<MT> normA(numev);

        // Back-transform the eigenvalues
        for (int i=0; i<numev; ++i) {
          MT mag2 = lapack.LAPY2(evals[i].realpart,evals[i].imagpart);
          mag2 = mag2*mag2;
          evals[i].realpart /=   mag2;
          evals[i].imagpart /= (-mag2);
        }

        // The problem is non-Hermitian.
        vector<int> curind(1);
        vector<MT> resnorm(1), tempnrm(1);
        Teuchos::RCP<const MV> Av_r, Av_i, Bv_r, Bv_i;
        Epetra_MultiVector Aevec(*Map,numev), Bevec(*Map,numev), res(*Map,1);

        // Compute A*evecs, B*evecs
        OPT::Apply( *A, *evecs, Aevec );
        OPT::Apply( *B, *evecs, Bevec );

        for (int i=0; i<numev;) {
          if (index[i]==0) {
            // Get views of the real part of A*evec,B*evec
            curind[0] = i;
            Av_r = MVT::CloneView( Aevec, curind );
            Bv_r = MVT::CloneView( Bevec, curind );

            // Compute set res = lambda*B*evec - A*evec
            // eigenvalue and eigenvector are both real
            MVT::MvAddMv(evals[i].realpart,*Bv_r,-1.0,*Av_r,res);

            // Compute the norm of the residual and increment counter
            MVT::MvNorm( res, resnorm );
            // Scale the residual
            normA[i] = resnorm[0];
            MT mag = MTT::magnitude(evals[i].realpart);
            if ( mag > MTT::one() ) {

//.........这里部分代码省略.........
        Mevecs = Teuchos::rcp( new Epetra_MultiVector(*Map,numev) );
        OPT::Apply( *M, *evecs, *Mevecs );
      }
      else {
        Mevecs = evecs;
      }
      
      Teuchos::SerialDenseMatrix<int,double> Breal(1,1), Bimag(1,1);
      while (i<numev) {
        if (index[i]==0) {
          // Get a view of the M*evecr
          curind[0] = i;
          tempMevec = MVT::CloneView( *Mevecs, curind );

          // Get a copy of A*evecr
          tempKevec = MVT::CloneCopy( Kevecs, curind );

          // Compute K*evecr - lambda*M*evecr
          Breal(0,0) = evals[i].realpart;
          MVT::MvTimesMatAddMv( -1.0, *tempMevec, Breal, 1.0, *tempKevec );

          // Compute the norm of the residual and increment counter
          MVT::MvNorm( *tempKevec, resnorm );
          normR[i] = resnorm[0]/Teuchos::ScalarTraits<double>::magnitude( evals[i].realpart );
          i++;
        } else {
          // Get a view of the real part of M*evecr
          curind[0] = i;
          tempMevec = MVT::CloneView( *Mevecs, curind );

          // Get a copy of K*evecr
          tempKevec = MVT::CloneCopy( Kevecs, curind );

          // Get a view of the imaginary part of the eigenvector (eveci)
          curind[0] = i+1;
          tempeveci = MVT::CloneView( *Mevecs, curind );

          // Set the eigenvalue into Breal and Bimag
          Breal(0,0) = evals[i].realpart;
          Bimag(0,0) = evals[i].imagpart;

          // Compute K*evecr - M*evecr*lambdar + M*eveci*lambdai
          MVT::MvTimesMatAddMv( -1.0, *tempMevec, Breal, 1.0, *tempKevec );
          MVT::MvTimesMatAddMv( 1.0, *tempeveci, Bimag, 1.0, *tempKevec );
          MVT::MvNorm( *tempKevec, tempnrm );

          // Get a copy of K*eveci
          tempKevec = MVT::CloneCopy( Kevecs, curind );

          // Compute K*eveci - M*eveci*lambdar - M*evecr*lambdai
          MVT::MvTimesMatAddMv( -1.0, *tempMevec, Bimag, 1.0, *tempKevec );
          MVT::MvTimesMatAddMv( -1.0, *tempeveci, Breal, 1.0, *tempKevec );
          MVT::MvNorm( *tempKevec, resnorm );

          // Compute the norms and scale by magnitude of eigenvalue
          normR[i] = lapack.LAPY2( tempnrm[i], resnorm[i] ) /
            lapack.LAPY2( evals[i].realpart, evals[i].imagpart );
          normR[i+1] = normR[i];

          i=i+2;
        }
      }
    }

    // Output computed eigenvalues and their direct residuals
    
    if (verbose && MyPID==0) {      
      cout.setf(std::ios_base::right, std::ios_base::adjustfield);	
      cout<<endl<< "Actual Residuals"<<endl;
      if (MyProblem->isHermitian()) {
        cout<< std::setw(16) << "Real Part"
          << std::setw(20) << "Direct Residual"<< endl;
        cout<<"-----------------------------------------------------------"<<endl;
        for (int i=0; i<numev; i++) {
          cout<< std::setw(16) << evals[i].realpart 
            << std::setw(20) << normR[i] << endl;
        }  
        cout<<"-----------------------------------------------------------"<<endl;
      } 
      else {
        cout<< std::setw(16) << "Real Part"
          << std::setw(16) << "Imag Part"
          << std::setw(20) << "Direct Residual"<< endl;
        cout<<"-----------------------------------------------------------"<<endl;
        for (int i=0; i<numev; i++) {
          cout<< std::setw(16) << evals[i].realpart 
            << std::setw(16) << evals[i].imagpart 
            << std::setw(20) << normR[i] << endl;
        }  
        cout<<"-----------------------------------------------------------"<<endl;
      }  
    }
  }

#ifdef EPETRA_MPI
    MPI_Finalize() ;
#endif
    return 0;

} // end BlockKrylovSchurEpetraExFile.cpp

//.........这里部分代码省略.........
    // Compute residuals.
    Teuchos::LAPACK<int,double> lapack;
    std::vector<double> normR(numev);

    // The problem is non-Hermitian.
    int i=0;
    std::vector<int> curind(1);
    std::vector<double> resnorm(1), tempnrm(1);
    Teuchos::RCP<MV> tempKevec, Mevecs;
    Teuchos::RCP<const MV> tempeveci, tempMevec;
    Epetra_MultiVector Kevecs(*Map,numev);

    // Compute K*evecs
    OPT::Apply( *K, *evecs, Kevecs );
    if (haveM) {
      Mevecs = Teuchos::rcp( new Epetra_MultiVector(*Map,numev) );
      OPT::Apply( *M, *evecs, *Mevecs );
    }
    else {
      Mevecs = evecs;
    }

    Teuchos::SerialDenseMatrix<int,double> Breal(1,1), Bimag(1,1);
    while (i<numev) {
      if (index[i]==0) {
        // Get a view of the M*evecr
        curind[0] = i;
        tempMevec = MVT::CloneView( *Mevecs, curind );

        // Get a copy of A*evecr
        tempKevec = MVT::CloneCopy( Kevecs, curind );

        // Compute K*evecr - lambda*M*evecr
        Breal(0,0) = evals[i].realpart;
        MVT::MvTimesMatAddMv( -1.0, *tempMevec, Breal, 1.0, *tempKevec );

        // Compute the norm of the residual and increment counter
        MVT::MvNorm( *tempKevec, resnorm );
        normR[i] = resnorm[0];
        i++;
      } else {
        // Get a view of the real part of M*evecr
        curind[0] = i;
        tempMevec = MVT::CloneView( *Mevecs, curind );

        // Get a copy of K*evecr
        tempKevec = MVT::CloneCopy( Kevecs, curind );

        // Get a view of the imaginary part of the eigenvector (eveci)
        curind[0] = i+1;
        tempeveci = MVT::CloneView( *Mevecs, curind );

        // Set the eigenvalue into Breal and Bimag
        Breal(0,0) = evals[i].realpart;
        Bimag(0,0) = evals[i].imagpart;

        // Compute K*evecr - M*evecr*lambdar + M*eveci*lambdai
        MVT::MvTimesMatAddMv( -1.0, *tempMevec, Breal, 1.0, *tempKevec );
        MVT::MvTimesMatAddMv( 1.0, *tempeveci, Bimag, 1.0, *tempKevec );
        MVT::MvNorm( *tempKevec, tempnrm );

        // Get a copy of K*eveci
        tempKevec = MVT::CloneCopy( Kevecs, curind );

        // Compute K*eveci - M*eveci*lambdar - M*evecr*lambdai
        MVT::MvTimesMatAddMv( -1.0, *tempMevec, Bimag, 1.0, *tempKevec );
        MVT::MvTimesMatAddMv( -1.0, *tempeveci, Breal, 1.0, *tempKevec );
        MVT::MvNorm( *tempKevec, resnorm );

        // Compute the norms and scale by magnitude of eigenvalue
        normR[i] = lapack.LAPY2( tempnrm[0], resnorm[0] );
        normR[i+1] = normR[i];

        i=i+2;
      }
    }

    // Output computed eigenvalues and their direct residuals
    if (verbose && MyPID==0) {
      cout.setf(std::ios_base::right, std::ios_base::adjustfield);
      cout<<endl<< "Actual Residuals"<<endl;
      cout<< std::setw(16) << "Real Part"
        << std::setw(16) << "Imag Part"
        << std::setw(20) << "Direct Residual"<< endl;
      cout<<"-----------------------------------------------------------"<<endl;
      for (int j=0; j<numev; j++) {
        cout<< std::setw(16) << evals[j].realpart
          << std::setw(16) << evals[j].imagpart
          << std::setw(20) << normR[j] << endl;
      }
      cout<<"-----------------------------------------------------------"<<endl;
    }
  }

#ifdef EPETRA_MPI
    MPI_Finalize() ;
#endif
    return 0;

} // end BlockKrylovSchurEpetraExFile.cpp