upload v4.3.0

Author: yoshiya-usui

doc/Manual_Of_FEMTIC.pdf

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

doc/Manual_Of_FEMTIC_v4.1.pdf

@@ -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 @@ class AnalysisControl{
 		// 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/AnalysisControl.cpp

@@ -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/AnalysisControl.h

@@ -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/CommonParameters.h

@@ -68,6 +68,10 @@ class DoubleSparseSquareMatrix : public DoubleSparseMatrix{
 	//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/DoubleSparseSquareMatrix.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