diff --git a/doc/Manual_Of_FEMTIC.pdf b/doc/Manual_Of_FEMTIC.pdf new file mode 100644 index 0000000..ebb82f4 Binary files /dev/null and b/doc/Manual_Of_FEMTIC.pdf differ diff --git a/doc/Manual_Of_FEMTIC_v4.1.pdf b/doc/Manual_Of_FEMTIC_v4.1.pdf deleted file mode 100644 index a2e866c..0000000 Binary files a/doc/Manual_Of_FEMTIC_v4.1.pdf and /dev/null differ diff --git a/src/AnalysisControl.cpp b/src/AnalysisControl.cpp index 9a492e0..6b4dcb2 100644 --- a/src/AnalysisControl.cpp +++ b/src/AnalysisControl.cpp @@ -1004,6 +1004,8 @@ void AnalysisControl::inputControlData(){ m_maxIterationIRWLSForLpOptimization = ibuf; inFile >> dbuf; m_thresholdIRWLSForLpOptimization = dbuf; + }else if (line.substr(0,19).compare("SENSE_MAT_DIRECTORY") == 0) { + inFile >> m_directoryOfOutOfCoreFilesForSensitivityMatrix; }else if( line.substr(0,3).compare("END") == 0 ){ break; }else{ @@ -1099,10 +1101,10 @@ void AnalysisControl::inputControlData(){ #ifdef _DEBUG_WRITE std::cout << "strEnv " << strEnv.str() << std::endl;// For debug #endif - if( putenv( strEnv.str().c_str() ) != 0 ){ - OutputFiles::m_logFile << "Error : Environment variable MKL_PARDISO_OOC_MAX_CORE_SIZE was not set correctly ! " << std::endl; - exit(1); - } + //if( putenv( strEnv.str().c_str() ) != 0 ){ + // OutputFiles::m_logFile << "Error : Environment variable MKL_PARDISO_OOC_MAX_CORE_SIZE was not set correctly ! " << std::endl; + // exit(1); + //} #endif OutputFiles::m_logFile << "# Maximum value of the memory used by out-of-core mode of forward solver : " << m_maxMemoryPARDISO << " [MB]" << std::endl; @@ -1444,6 +1446,10 @@ void AnalysisControl::inputControlData(){ OutputFiles::m_logFile << "# Maximum number of retrials : " << m_numCutbackMax << "." << std::endl; } + if (!getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { + OutputFiles::m_logFile << "# Directory of out-of-core files for the sensitivitry matrix: " << getDirectoryOfOutOfCoreFilesForSensitivityMatrix() << std::endl; + } + OutputFiles::m_logFile << "#==============================================================================" << std::endl; } @@ -1750,6 +1756,11 @@ double AnalysisControl::getThresholdIRWLSForLpOptimization() const{ return m_thresholdIRWLSForLpOptimization; } +// Get directory of out-of-core files for the sensitivitry matrix +std::string AnalysisControl::getDirectoryOfOutOfCoreFilesForSensitivityMatrix() const { + return m_directoryOfOutOfCoreFilesForSensitivityMatrix; +} + #ifdef _ANISOTOROPY // Get type of anisotropy int AnalysisControl::getTypeOfAnisotropy() const{ diff --git a/src/AnalysisControl.h b/src/AnalysisControl.h index 4693e2e..ac68d73 100644 --- a/src/AnalysisControl.h +++ b/src/AnalysisControl.h @@ -273,6 +273,9 @@ class AnalysisControl{ // Get threshold value for deciding convergence about IRWLS for Lp optimization double getThresholdIRWLSForLpOptimization() const; + // Get directory of out-of-core files for the sensitivitry matrix + std::string getDirectoryOfOutOfCoreFilesForSensitivityMatrix() const; + #ifdef _ANISOTOROPY // Get type of anisotropy int getTypeOfAnisotropy() const; @@ -559,6 +562,9 @@ private: // Threshold value for deciding convergence about IRWLS for Lp optimization double m_thresholdIRWLSForLpOptimization; + // Directory of out-of-core files for the sensitivitry matrix + std::string m_directoryOfOutOfCoreFilesForSensitivityMatrix; + #ifdef _ANISOTOROPY // Type of anisotropy int m_typeOfAnisotropy; diff --git a/src/CommonParameters.h b/src/CommonParameters.h index 9cb68ed..fd39f72 100644 --- a/src/CommonParameters.h +++ b/src/CommonParameters.h @@ -162,7 +162,7 @@ static char programName[]="femtic"; // [MajorVersion#].[MinorVersion#].[Revision#] // x.x.xa -> alpha version // x.x.xb -> beta version -static char versionID[]="4.2.5"; +static char versionID[] = "4.3.0"; } diff --git a/src/DoubleSparseSquareMatrix.cpp b/src/DoubleSparseSquareMatrix.cpp index 163f2ae..2994486 100644 --- a/src/DoubleSparseSquareMatrix.cpp +++ b/src/DoubleSparseSquareMatrix.cpp @@ -27,6 +27,14 @@ #include "DoubleSparseSquareMatrix.h" #include "OutputFiles.h" #include +#include + +#ifdef _DEBUG_WRITE_FOR_BOTTOM_RESISTIVITY +#ifdef _LINUX +#include +#include +#endif +#endif //Default Constructer DoubleSparseSquareMatrix::DoubleSparseSquareMatrix(): @@ -135,6 +143,147 @@ void DoubleSparseSquareMatrix::solvePhaseMatrixSolver( const int nrhs, double* r m_pardisoSolver.solve( m_rowIndex, m_columns, m_values, nrhs, rhs, solution ); } +//Solve phase of matrix solver by the conjugate gradient method with the point Jacobi preconditioner +//@note Matrix should be symmetric +void DoubleSparseSquareMatrix::solvePhaseMatrixSolverByPCGPointJacobi(const int nrhs, double* rhs, double* solution) const{ + assert(m_hasConvertedToCRSFormat); + + const int maxIterationNumber = m_numRows; + const double eps = 1.0e-20; + double* invDiagonals = new double[m_numRows]; + double* workP = new double[m_numRows]; + double* workR = new double[m_numRows];// Residuals + double* workQ = new double[m_numRows]; + double* workX = new double[m_numRows];// Solution vector + double* workZ = new double[m_numRows]; + + for (int irow = 0; irow < m_numRows; ++irow) + { + for (int j = m_rowIndex[irow]; j < m_rowIndex[irow + 1]; ++j) + { + if (irow == m_columns[j]) + { + invDiagonals[irow] = 1.0 / m_values[j]; + } + } + } + for (int irhs = 0; irhs < nrhs; ++irhs) + { + // Initial solution is a zero vector + for (int irow = 0; irow < m_numRows; ++irow) + { + workX[irow] = 0.0; + } + // [r0] = [b] - [A][x0] + double normOfRhsVector(0.0); + for (int irow = 0; irow < m_numRows; ++irow) + { + const long long int index = static_cast(irow) + static_cast(irhs) * static_cast(m_numRows); + normOfRhsVector += rhs[index] * rhs[index]; + workR[irow] = rhs[index]; + } + int iter = 0; + double rhoPre(0.0); + for (; iter < maxIterationNumber; ++iter) + { + // [z] = [M]^-1[r] + for (int irow = 0; irow < m_numRows; ++irow) + { + workZ[irow] = invDiagonals[irow] * workR[irow]; + } + // rho = [r]T[z] + double rho(0.0); + for (int irow = 0; irow < m_numRows; ++irow) + { + rho += workR[irow] * workZ[irow]; + } + if (iter == 0) + { + // [p0] - [z0] + for (int irow = 0; irow < m_numRows; ++irow) + { + workP[irow] = workZ[irow]; + } + } + else + { + // [p] = [z] + beta*[p] + const double beta = rho / rhoPre; + for (int irow = 0; irow < m_numRows; ++irow) + { + workP[irow] = workZ[irow] + beta * workP[irow]; + } + } + // [q] = [A][p] + for (int irow = 0; irow < m_numRows; ++irow) + { + workQ[irow] = 0.0; + for (int j = m_rowIndex[irow]; j < m_rowIndex[irow + 1]; ++j) + { + workQ[irow] += m_values[j] * workP[m_columns[j]]; + } + } + // alpha = rho / [p]T[q] + double pq(0.0); + for (int irow = 0; irow < m_numRows; ++irow) + { + pq += workP[irow] * workQ[irow]; + } + const double alpha = rho / pq; + // [x] = [x] + alpha * [p] + // [r] = [r] - alpha * [q] + for (int irow = 0; irow < m_numRows; ++irow) + { + workX[irow] += alpha * workP[irow]; + workR[irow] -= alpha * workQ[irow]; + } + // Check convergence + double normOfResidualVector(0.0); + for (int irow = 0; irow < m_numRows; ++irow) + { + normOfResidualVector += workR[irow] * workR[irow]; + } + if( sqrt(normOfResidualVector/ normOfRhsVector) < eps ) + { + break; + } + rhoPre = rho; + } + if (iter >= maxIterationNumber) { + OutputFiles::m_logFile << "Error : PCG solver is not converged !!" << std::endl; + exit(1); + } + else { + OutputFiles::m_logFile << "# PCG solver is converged after " << iter << " iterations." << std::endl; + } + for (int irow = 0; irow < m_numRows; ++irow) + { + const long long int index = static_cast(irow) + static_cast(irhs) * static_cast(m_numRows); + solution[index] = workX[irow]; + } + } + +#ifdef _DEBUG_WRITE_FOR_BOTTOM_RESISTIVITY +#ifdef _LINUX + { + struct rusage r; + if (getrusage(RUSAGE_SELF, &r) != 0) { + /*Failure*/ + } + OutputFiles::m_logFile << "maxrss= " << r.ru_maxrss << std::endl; + } +#endif +#endif + + delete[] invDiagonals; + delete[] workP; + delete[] workR; + delete[] workQ; + delete[] workX; + delete[] workZ; + +} + //Release memory of matrix solver void DoubleSparseSquareMatrix::releaseMemoryMatrixSolver(){ if( m_pardisoSolver.getSolutionStage() > PARDISOSolver::MEMORY_RELEASED ){ diff --git a/src/DoubleSparseSquareMatrix.h b/src/DoubleSparseSquareMatrix.h index cf81e6a..7c7b619 100644 --- a/src/DoubleSparseSquareMatrix.h +++ b/src/DoubleSparseSquareMatrix.h @@ -68,6 +68,10 @@ public: //Solve phase of matrix solver void solvePhaseMatrixSolver( const int nrhs, double* rhs, double* solution ); + //Solve phase of matrix solver by the conjugate gradient method with the point Jacobi preconditioner + //@note Matrix should be symmetric + void solvePhaseMatrixSolverByPCGPointJacobi(const int nrhs, double* rhs, double* solution) const; + //Release memory of matrix solver void releaseMemoryMatrixSolver(); diff --git a/src/Inversion.cpp b/src/Inversion.cpp index 2332a52..722b662 100644 --- a/src/Inversion.cpp +++ b/src/Inversion.cpp @@ -185,6 +185,13 @@ void Inversion::calculateSensitivityMatrix( const int freqIDAmongThisPE, const d // Write sensitivity matrix to out-of-core file --- //------------------------------------------------- std::ostringstream fileName; + if ( !ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty() ){ +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqIDGlobal; //std::ofstream outputFile(fileName.str().c_str(), std::ios::out|std::ios::binary|std::ios::trunc); //if( outputFile.fail() ){ @@ -713,10 +720,18 @@ void Inversion::multiplyModelTransformingJacobian( const int numData, const int void Inversion::deleteOutOfCoreFileAll(){ const ObservedData* const ptrObservedData = ObservedData::getInstance(); + const AnalysisControl* const ptrAnalysisControl = AnalysisControl::getInstance(); const int nFreq = ptrObservedData->getNumOfFrequenciesCalculatedByThisPE(); for( int iFreq = 0; iFreq < nFreq; ++iFreq ){ const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(iFreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; FILE* fp = fopen( fileName.str().c_str(), "rb" ); if( fp != NULL ){// File exists diff --git a/src/InversionGaussNewtonDataSpace.cpp b/src/InversionGaussNewtonDataSpace.cpp index 077b76e..fb645f9 100644 --- a/src/InversionGaussNewtonDataSpace.cpp +++ b/src/InversionGaussNewtonDataSpace.cpp @@ -32,6 +32,13 @@ #include #include +#ifdef _DEBUG_WRITE_FOR_BOTTOM_RESISTIVITY +#ifdef _LINUX +#include +#include +#endif +#endif + #include "mkl_cblas.h" #include "mkl_lapacke.h" #include "mpi.h" @@ -71,7 +78,7 @@ void InversionGaussNewtonDataSpace::inversionCalculation(){ } // Perform inversion by the new method -void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ +void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const { const bool useBLAS = true; @@ -84,7 +91,7 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ std::ostringstream oss; oss << "debug_" << myProcessID << ".txt"; std::ofstream fout; - fout.open( oss.str().c_str() ); + fout.open(oss.str().c_str()); #endif ObservedData* const ptrObservedData = ObservedData::getInstance(); @@ -100,11 +107,22 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ // Calculate constraining matrix //--------------------------------- OutputFiles::m_logFile << "# Calculate constraining matrix. " << ptrAnalysisControl->outputElapsedTime() << std::endl; +#ifdef _DEBUG_WRITE_FOR_BOTTOM_RESISTIVITY +#ifdef _LINUX + { + struct rusage r; + if (getrusage(RUSAGE_SELF, &r) != 0) { + /*Failure*/ + } + OutputFiles::m_logFile << "maxrss= " << r.ru_maxrss << std::endl; +} +#endif +#endif RougheningSquareMatrix constrainingMatrix; - calcConstrainingMatrix( constrainingMatrix ); + calcConstrainingMatrix(constrainingMatrix); RougheningSquareMatrix transposedConstrainingMatrix; constrainingMatrix.makeTansposedMatrix( transposedConstrainingMatrix ); -//#ifdef _DEBUG_WRITE +//#ifdef _DEBUG_WRITE_FOR_BOTTOM_RESISTIVITY // constrainingMatrix.calcSingularValues(); //#endif @@ -127,14 +145,49 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ constrainingMatrix.calcMatrixVectorProductUsingTransposedMatrix( vectorRxm, dVector ); } +#ifdef _PCG +#ifdef _PCG + const bool usePCG = true; +#endif + if (!usePCG) { + //---------------------------------- + // Initialization + //---------------------------------- + std::ostringstream oocHeaderName; + oocHeaderName << "ooc_temp_inv_3D_PE" << myProcessID; + constrainingMatrix.initializeMatrixSolver(oocHeaderName.str(), ptrAnalysisControl->getModeOfPARDISO()); + oocHeaderName << "_trans"; + transposedConstrainingMatrix.initializeMatrixSolver(oocHeaderName.str(), ptrAnalysisControl->getModeOfPARDISO()); + + //---------------------------------- + // Analysis + //---------------------------------- + OutputFiles::m_logFile << "# Analyse constraining matrix." << ptrAnalysisControl->outputElapsedTime() << std::endl; + constrainingMatrix.analysisPhaseMatrixSolver(); + OutputFiles::m_logFile << "# Analyse transposed constraining matrix." << ptrAnalysisControl->outputElapsedTime() << std::endl; + transposedConstrainingMatrix.analysisPhaseMatrixSolver(); + + //---------------------------------- + // Factorization + //---------------------------------- + OutputFiles::m_logFile << "# Factorize constraining matrix." << ptrAnalysisControl->outputElapsedTime() << std::endl; + constrainingMatrix.factorizationPhaseMatrixSolver(); + OutputFiles::m_logFile << "# Factorize transposed constraining matrix." << ptrAnalysisControl->outputElapsedTime() << std::endl; + transposedConstrainingMatrix.factorizationPhaseMatrixSolver(); + } + else + { + OutputFiles::m_logFile << "# PCG solver is used for inverting the roughening matrix." << ptrAnalysisControl->outputElapsedTime() << std::endl; + } +#else //---------------------------------- // Initialization //---------------------------------- std::ostringstream oocHeaderName; oocHeaderName << "ooc_temp_inv_3D_PE" << myProcessID; - constrainingMatrix.initializeMatrixSolver( oocHeaderName.str(), ptrAnalysisControl->getModeOfPARDISO() ); + constrainingMatrix.initializeMatrixSolver(oocHeaderName.str(), ptrAnalysisControl->getModeOfPARDISO()); oocHeaderName << "_trans"; - transposedConstrainingMatrix.initializeMatrixSolver( oocHeaderName.str(), ptrAnalysisControl->getModeOfPARDISO() ); + transposedConstrainingMatrix.initializeMatrixSolver(oocHeaderName.str(), ptrAnalysisControl->getModeOfPARDISO()); //---------------------------------- // Analysis @@ -151,6 +204,18 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ constrainingMatrix.factorizationPhaseMatrixSolver(); OutputFiles::m_logFile << "# Factorize transposed constraining matrix." << ptrAnalysisControl->outputElapsedTime() << std::endl; transposedConstrainingMatrix.factorizationPhaseMatrixSolver(); +#endif +#ifdef _DEBUG_WRITE_FOR_BOTTOM_RESISTIVITY +#ifdef _LINUX + { + struct rusage r; + if (getrusage(RUSAGE_SELF, &r) != 0) { + /*Failure*/ + } + OutputFiles::m_logFile << "maxrss= " << r.ru_maxrss << std::endl; + } +#endif +#endif //------------------------------------------------------------------------------------------------------ // Calculate residual vector and sensitivity matrix multiplied by inverse of the constraing matrix @@ -182,6 +247,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(iFreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileName.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; @@ -259,9 +331,24 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ const int numRHSDivided = iDiv < numAdds ? numRHSDividedWithoutOdds + 1 : numRHSDividedWithoutOdds; OutputFiles::m_logFile << "# Solve phase is performed simultaneously for " << numRHSDivided << " right-hand sides" << ptrAnalysisControl->outputElapsedTime() << std::endl; const long long index = static_cast(numModel) * iRhsStart; +#ifdef _PCG + if (usePCG) { + transposedConstrainingMatrix.solvePhaseMatrixSolverByPCGPointJacobi(numRHSDivided, &sensitivityMatrix[index], &sensitivityMatrixTemp[index]);// Solve + } + else { + transposedConstrainingMatrix.solvePhaseMatrixSolver(numRHSDivided, &sensitivityMatrix[index], &sensitivityMatrixTemp[index]);// Solve + } +#else transposedConstrainingMatrix.solvePhaseMatrixSolver( numRHSDivided, &sensitivityMatrix[index], &sensitivityMatrixTemp[index] );// Solve +#endif iRhsStart += numRHSDivided; } +#ifdef _DEBUG_WRITE + for (long long int index = 0; index < numComps; ++index) + { + fout << "sensitivityMatrixTemp " << index << " " << sensitivityMatrixTemp[index] << std::endl; + } +#endif //------------------------------------------------------------------------------------ // Write sensitivity matrix multiplied by inverse of constraining matrix @@ -330,7 +417,23 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ //------------------------------------------ double* vectorInvRTxJTxResidial = new double[numModel]; OutputFiles::m_logFile << "# Solve phase." << ptrAnalysisControl->outputElapsedTime() << std::endl; +#ifdef _PCG + if (usePCG) { + transposedConstrainingMatrix.solvePhaseMatrixSolverByPCGPointJacobi(1, vectorJTxResidialGlobal, vectorInvRTxJTxResidial); + } + else + { + transposedConstrainingMatrix.solvePhaseMatrixSolver(1, vectorJTxResidialGlobal, vectorInvRTxJTxResidial); + } +#ifdef _DEBUG_WRITE + for (int index = 0; index < numModel; ++index) + { + fout << "vectorInvRTxJTxResidial " << index << " " << vectorInvRTxJTxResidial[index] << std::endl; + } +#endif +#else transposedConstrainingMatrix.solvePhaseMatrixSolver( 1, vectorJTxResidialGlobal, vectorInvRTxJTxResidial ); +#endif delete [] vectorJTxResidialGlobal; #ifdef _DEBUG_WRITE for( int iMdl = 0; iMdl < numModel; ++iMdl ){ @@ -357,7 +460,23 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ // vectorRxm : Rm + inv[R]T*[J]T*rd // vectorInvRTRd : inv[R]*( Rm + inv[R]T*[J]T*rd ) //------------------------------------------------- +#ifdef _PCG + if (usePCG) { + constrainingMatrix.solvePhaseMatrixSolverByPCGPointJacobi(1, vectorRxm, vectorInvRTRd); + } + else + { + constrainingMatrix.solvePhaseMatrixSolver(1, vectorRxm, vectorInvRTRd); + } +#ifdef _DEBUG_WRITE + for (int index = 0; index < numModel; ++index) + { + fout << "vectorInvRTRd " << index << " " << vectorInvRTRd[index] << std::endl; + } +#endif +#else constrainingMatrix.solvePhaseMatrixSolver( 1, vectorRxm, vectorInvRTRd ); +#endif delete [] vectorRxm; #ifdef _DEBUG_WRITE @@ -383,6 +502,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(iFreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileName.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; @@ -521,6 +647,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ int offsetRows(0); for( int ifreqLeft = 0 ;ifreqLeft < nFreq; ++ifreqLeft ){ std::ostringstream fileNameLeft; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileNameLeft << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileNameLeft << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileNameLeft << "sensMatFreq" << ifreqLeft << "Mod"; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileNameLeft.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; @@ -536,6 +669,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ int offsetCols = offsetRows; for( int ifreqRight = ifreqLeft ;ifreqRight < nFreq; ++ifreqRight ){ std::ostringstream fileNameRight; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileNameRight << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileNameRight << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileNameRight << "sensMatFreq" << ifreqRight << "Mod"; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileNameRight.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; @@ -742,6 +882,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ assert( offset == ptrObservedData->getNumObservedDataThisPEAccumulated( iFreq ) ); const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(iFreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileName.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; int numDataThisFreq(0); @@ -838,10 +985,30 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ // resultVector : inv[ [R]T*[R] ] * [ RTRm + [J]T*rd ] - inv[ [R]T*[R] ] * [J]T * inv[ [I] + [J] * inv[ [R]T*[R] ] *[J]T ] * [J] * inv[ [R]T*[R] ] * [ RTRm + [J]T*rd ] //------------------------------------------------------------------------------------------------------------------- double* tempVector = new double[numModel]; +#ifdef _PCG + if (usePCG) { + transposedConstrainingMatrix.solvePhaseMatrixSolverByPCGPointJacobi(1, dVector, tempVector); + constrainingMatrix.solvePhaseMatrixSolverByPCGPointJacobi(1, tempVector, dVector); + } + else + { + transposedConstrainingMatrix.solvePhaseMatrixSolver(1, dVector, tempVector); + constrainingMatrix.solvePhaseMatrixSolver(1, tempVector, dVector); + } +#ifdef _DEBUG_WRITE + for (int index = 0; index < numModel; ++index) + { + fout << "tempVector " << index << " " << tempVector[index] << std::endl; + } +#endif + transposedConstrainingMatrix.releaseMemory(); + constrainingMatrix.releaseMemory(); +#else transposedConstrainingMatrix.solvePhaseMatrixSolver( 1, dVector, tempVector ); transposedConstrainingMatrix.releaseMemory(); constrainingMatrix.solvePhaseMatrixSolver( 1, tempVector, dVector ); constrainingMatrix.releaseMemory(); +#endif delete [] tempVector; //------------------------------------------------------------------------------------------------------------------- } @@ -878,6 +1045,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ for( int ifreq = 0 ;ifreq < nFreqThisPE; ++ifreq ){ const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(ifreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; fileName << "Mod"; if( remove( fileName.str().c_str() ) != 0 ){ @@ -1015,6 +1189,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(iFreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileName.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; @@ -1071,7 +1252,7 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa //------------------------------------------------------------------------------------- // Calculate sensitivity matrix multiplied by inverse of constraining matrix //------------------------------------------------------------------------------------- - // sensitivityMatrix : JT => inv[R]*[J]T + // sensitivityMatrix : JT => inv([R]T[R])*[J]T //------------------------------------------------------------------------------------- OutputFiles::m_logFile << "# Multiply transposed sensitivity matrix by inverse of constraining matrix. " << ptrAnalysisControl->outputElapsedTime() << std::endl; const long long numComps = static_cast(numDataThisFreq) * static_cast(numModel); @@ -1192,6 +1373,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(iFreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileName.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; @@ -1330,6 +1518,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa int offsetRows(0); for( int ifreqLeft = 0 ;ifreqLeft < nFreq; ++ifreqLeft ){ std::ostringstream fileNameLeft; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileNameLeft << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileNameLeft << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileNameLeft << "sensMatFreq" << ifreqLeft; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileNameLeft.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; @@ -1345,6 +1540,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa int offsetCols = offsetRows; for( int ifreqRight = ifreqLeft ;ifreqRight < nFreq; ++ifreqRight ){ std::ostringstream fileNameRight; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileNameRight << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileNameRight << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileNameRight << "sensMatFreq" << ifreqRight << "Mod"; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileNameRight.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; @@ -1551,6 +1753,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa assert( offset == ptrObservedData->getNumObservedDataThisPEAccumulated( iFreq ) ); const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(iFreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; OutputFiles::m_logFile << "# Read sensitivity matrix from "<< fileName.str() << "." << ptrAnalysisControl->outputElapsedTime() << std::endl; int numDataThisFreq(0); @@ -1686,6 +1895,13 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa for( int ifreq = 0 ;ifreq < nFreqThisPE; ++ifreq ){ const int freqID = ptrObservedData->getIDsOfFrequenciesCalculatedByThisPE(ifreq); std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; fileName << "Mod"; if( remove( fileName.str().c_str() ) != 0 ){ diff --git a/src/InversionGaussNewtonModelSpace.cpp b/src/InversionGaussNewtonModelSpace.cpp index e98f31c..aa97b9e 100644 --- a/src/InversionGaussNewtonModelSpace.cpp +++ b/src/InversionGaussNewtonModelSpace.cpp @@ -167,6 +167,13 @@ void InversionGaussNewtonModelSpace::inversionCalculation(){ const int freqID = iFreq; std::ostringstream fileName; + if (!ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix().empty()) { +#ifdef _LINUX + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\/"; +#else + fileName << ptrAnalysisControl->getDirectoryOfOutOfCoreFilesForSensitivityMatrix() + "\\"; +#endif + } fileName << "sensMatFreq" << freqID; FILE* fp = fopen( fileName.str().c_str(), "rb" ); if( fp == NULL ){