Update README.md

README.md

@@ -7,12 +7,12 @@ FEMTIC is a 3-D magnetotelluric inversion code based on the following studies. F
 
 *Y. Usui, T. Kasaya, Y. Ogawa and H. Iwamoto, Marine magnetotelluric inversion with an unstructured tetrahedral mesh, Geophys. J. Int., 214(2): 952-974, https://doi.org/10.1093/gji/ggy171, 2018.*
 
-*Y. Usui, Applicability evaluation of non-conforming deformed hexahedral mesh for marine magnetotellurics, Japan Geoscience Union Meeting 2021, 2021*
+*Y. Usui, M. Uyeshima, H. Hase, H. Ichihara, K. Aizawa, T. Koyama, et al. Three-dimensional electrical resistivity structure beneath a strain concentration area in the back-arc side of the northeastern Japan arc. Journal of Geophysical Research: Solid Earth, 129, e2023JB028522. http://doi.org/10.1029/2023JB028522, 2024.*
 
 **The website of FEMTIC:**
 https://sites.google.com/view/yoshiyausui/femtic
 
-# Functional overview
+## Functional overview
 FEMTIC gives a three-dimensional electrical resistivity structure from the response functions at observation points on the Earth's surface.
 
 **Mesh type**: Tetrahedral mesh / Hexahedral brick mesh / Non-conforming deformed hexahedral mesh
@@ -28,7 +28,11 @@ FEMTIC gives a three-dimensional electrical resistivity structure from the respo
 **Regularization**: L2 regularization with Laplacian filter / L2 regularization with difference filter /  L1 regularization with difference filter
 
 
-# Release note
+## Release note
+***v4.3*** Feb. 6, 2025: The directory for the out-of-core files of the sensitivity matrix becomes changeable.
+
+***v4.2*** Mar. 25, 2024: I modified some parts to allow the use of large-scale models and large datasets.
+
 ***v4.1*** Nov. 9, 2021: This new version supports difference filter. The error calculation of log10(apparent resistivity) is modified. Rotation angles of distortion matrix are limited to from -90 to 90 (deg.) when gains and rotations of galvanic distortion are estimated.
 
 ***v4.0*** Jun. 3, 2021: This new version supports non-conforming deformed hexahedral mesh.
@@ -39,7 +43,7 @@ FEMTIC gives a three-dimensional electrical resistivity structure from the respo
 
 ***v3.4.6*** Sep. 2, 2020: This version allows us to make resistivity of every individual subsurface element to be a different model parameter, in analogy with other 3-D inversion code.
 
-# Pre/post-processing tools for FEMTIC
+## Pre/post-processing tools for FEMTIC
 Some pre/post-processing tools for FEMTIC, including meshing tools, and their manuals can be downloaded from GitHub. Results of FEMTIC can be visualized by [ParaView](https://www.paraview.org/).
 
 [makeDHexaMesh](https://github.com/yoshiya-usui/makeDHexaMesh.git): Tool for making non-conforming deformed hexahedral mesh for FEMTIC

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{

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;

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";
 
 }
 

src/DoubleSparseSquareMatrix.cpp

@@ -27,6 +27,14 @@
 #include "DoubleSparseSquareMatrix.h"
 #include "OutputFiles.h"
 #include <assert.h>
+#include <math.h>
+
+#ifdef _DEBUG_WRITE_FOR_BOTTOM_RESISTIVITY
+#ifdef _LINUX
+#include <sys/time.h>
+#include <sys/resource.h>
+#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<long long int>(irow) + static_cast<long long int>(irhs) * static_cast<long long int>(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<long long int>(irow) + static_cast<long long int>(irhs) * static_cast<long long int>(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 ){

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();
 

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

src/InversionGaussNewtonDataSpace.cpp

@@ -32,6 +32,13 @@
 #include <string.h>
 #include <assert.h>
 
+#ifdef _DEBUG_WRITE_FOR_BOTTOM_RESISTIVITY
+#ifdef _LINUX
+#include <sys/time.h>
+#include <sys/resource.h>
+#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<long long>(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<long long>(numDataThisFreq) * static_cast<long long>(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 ){

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 ){