From 8eaec75b2e7d533bf5c855deab6e4c0f6a1ee5eb Mon Sep 17 00:00:00 2001 From: yoshiya-usui Date: Wed, 1 Mar 2023 06:12:56 +0900 Subject: [PATCH] modify some files for large-scale models --- src/ComplexSparseMatrix.cpp | 37 +++++++++++++++------------ src/ComplexSparseSquareMatrix.cpp | 4 +-- src/ComplexSparseSquareMatrix.h | 2 +- src/DoubleSparseSquareMatrix.cpp | 4 +-- src/DoubleSparseSquareMatrix.h | 2 +- src/Forward3D.cpp | 6 ++--- src/InversionGaussNewtonDataSpace.cpp | 32 +++++++++++------------ 7 files changed, 45 insertions(+), 42 deletions(-) diff --git a/src/ComplexSparseMatrix.cpp b/src/ComplexSparseMatrix.cpp index 99c7758..b5bf232 100644 --- a/src/ComplexSparseMatrix.cpp +++ b/src/ComplexSparseMatrix.cpp @@ -61,11 +61,11 @@ ComplexSparseMatrix::ComplexSparseMatrix( const int nrows, const int ncols, cons assert( ncols > 0 ); assert( nrhs > 0 ); - const int num = m_numRows*m_numRightHandSideVectors; + const long long num = static_cast(m_numRows)*static_cast(m_numRightHandSideVectors); m_rightHandSideVector = new std::complex[num]; - for( int i = 0; i < num; ++i ){ + for( long long i = 0; i < num; ++i ){ m_rightHandSideVector[i] = std::complex(0.0,0.0); // Initialize - } + } } // Destructer @@ -125,9 +125,9 @@ void ComplexSparseMatrix::setNumRowsAndColumns( const int nrows, const int ncols m_rightHandSideVector = NULL; } - const int num = m_numRows*m_numRightHandSideVectors; + const long long num = static_cast(m_numRows)*static_cast(m_numRightHandSideVectors); m_rightHandSideVector = new std::complex[num]; - for( int i = 0; i < num; ++i ){ + for( long long i = 0; i < num; ++i ){ m_rightHandSideVector[i] = std::complex(0.0,0.0); // Initialize } @@ -374,14 +374,16 @@ void ComplexSparseMatrix::addRightHandSideVector( const int row, const std::comp assert( irhs <= m_numRightHandSideVectors - 1 ); assert( irhs >= 0 ); - m_rightHandSideVector[ row + m_numRows * irhs ] += val; + const long long index = static_cast(row) + static_cast(m_numRows)*static_cast(irhs); + m_rightHandSideVector[index] += val; + } //Zero clear non-zero values of the right hand side vector void ComplexSparseMatrix::zeroClearRightHandSideVector(){ - const int num = m_numRows*m_numRightHandSideVectors; - for( int i = 0; i < num; ++i ){ + const long long num = static_cast(m_numRows)*static_cast(m_numRightHandSideVectors); + for( long long i = 0; i < num; ++i ){ m_rightHandSideVector[i] = std::complex(0.0,0.0); // Zero clear } @@ -444,9 +446,9 @@ void ComplexSparseMatrix::initializeMatrixAndRhsVectors( const int nrows, const //m_matrixTripletFormat = new std::set[nrows]; m_matrixTripletFormat = new std::map< int, std::complex >[nrows]; - const int num = m_numRows*m_numRightHandSideVectors; + const long long num = static_cast(m_numRows)*static_cast(m_numRightHandSideVectors); m_rightHandSideVector = new std::complex[num]; - for( int i = 0; i < num; ++i ){ + for( long long i = 0; i < num; ++i ){ m_rightHandSideVector[i] = std::complex(0.0,0.0); // Initialize } @@ -530,9 +532,9 @@ void ComplexSparseMatrix::reallocateMemoryForRightHandSideVectors( const int nrh m_numRightHandSideVectors = nrhs; - const int num = m_numRows*m_numRightHandSideVectors; + const long long num = static_cast(m_numRows)*static_cast(m_numRightHandSideVectors); m_rightHandSideVector = new std::complex[num]; - for( int i = 0; i < num; ++i ){ + for( long long i = 0; i < num; ++i ){ m_rightHandSideVector[i] = std::complex(0.0,0.0); // Initialize } @@ -574,17 +576,18 @@ void ComplexSparseMatrix::copyRhsVector( std::complex* vecOut ) const{ // } //} - memcpy( vecOut, m_rightHandSideVector, sizeof(std::complex)*(m_numRightHandSideVectors*m_numRows) ); + const long long num = static_cast(m_numRightHandSideVectors)*static_cast(m_numRows); + memcpy( vecOut, m_rightHandSideVector, static_cast(sizeof(std::complex))*num ); } //Copy specified components of right-hand side vector to another vector void ComplexSparseMatrix::copyRhsVector( const int numCompsCopied, const int* const compsCopied, std::complex* vecOut ) const{ - for( int j = 0; j < m_numRightHandSideVectors; ++j ){ - const int offset = m_numRows * j; - for( int i = 0; i < numCompsCopied; ++i ){ - vecOut[ i + offset ] = m_rightHandSideVector[ compsCopied[i] + offset ]; + for( long long j = 0; j < m_numRightHandSideVectors; ++j ){ + const long long offset = static_cast(m_numRows) * j; + for( long long i = 0; i < numCompsCopied; ++i ){ + vecOut[ i + offset ] = m_rightHandSideVector[ static_cast(compsCopied[i]) + offset ]; } } diff --git a/src/ComplexSparseSquareMatrix.cpp b/src/ComplexSparseSquareMatrix.cpp index 0cb5e74..d6592f0 100644 --- a/src/ComplexSparseSquareMatrix.cpp +++ b/src/ComplexSparseSquareMatrix.cpp @@ -333,11 +333,11 @@ void ComplexSparseSquareMatrix::factorizationPhaseMatrixSolver(){ } //Solve phase of matrix solver with a specified number of right-hand side -void ComplexSparseSquareMatrix::solvePhaseMatrixSolver( std::complex* solution, const int iRhsStart ,const int nRhs ){ +void ComplexSparseSquareMatrix::solvePhaseMatrixSolver( std::complex* solution, const long long iRhsStart ,const int nRhs ){ assert( m_hasConvertedToCRSFormat ); - const long long index = static_cast(m_numRows) * static_cast(iRhsStart); + const long long index = static_cast(m_numRows) * iRhsStart; m_pardisoSolver.solve( m_rowIndex, m_columns, m_values, nRhs, &m_rightHandSideVector[index], solution ); } diff --git a/src/ComplexSparseSquareMatrix.h b/src/ComplexSparseSquareMatrix.h index d26930c..13a7f3d 100644 --- a/src/ComplexSparseSquareMatrix.h +++ b/src/ComplexSparseSquareMatrix.h @@ -79,7 +79,7 @@ public: void factorizationPhaseMatrixSolver(); //Solve phase of matrix solver with a specified number of right-hand side - void solvePhaseMatrixSolver( std::complex* solution, const int iRhsStart ,const int nRhs ); + void solvePhaseMatrixSolver( std::complex* solution, const long long iRhsStart ,const int nRhs ); //Solve phase of matrix solver void solvePhaseMatrixSolver( std::complex* solution ); diff --git a/src/DoubleSparseSquareMatrix.cpp b/src/DoubleSparseSquareMatrix.cpp index e4c6dee..163f2ae 100644 --- a/src/DoubleSparseSquareMatrix.cpp +++ b/src/DoubleSparseSquareMatrix.cpp @@ -117,9 +117,9 @@ void DoubleSparseSquareMatrix::factorizationPhaseMatrixSolver(){ } //Solve phase of matrix solver with a specified number of right-hand side -void DoubleSparseSquareMatrix::solvePhaseMatrixSolver( double* solution, const int iRhsStart ,const int nRhs ){ +void DoubleSparseSquareMatrix::solvePhaseMatrixSolver( double* solution, const long long iRhsStart ,const int nRhs ){ assert( m_hasConvertedToCRSFormat ); - const long long index = static_cast(m_numRows) * static_cast(iRhsStart); + const long long index = static_cast(m_numRows) * iRhsStart; m_pardisoSolver.solve( m_rowIndex, m_columns, m_values, nRhs, &m_rightHandSideVector[index], solution ); } diff --git a/src/DoubleSparseSquareMatrix.h b/src/DoubleSparseSquareMatrix.h index 42c036d..cf81e6a 100644 --- a/src/DoubleSparseSquareMatrix.h +++ b/src/DoubleSparseSquareMatrix.h @@ -60,7 +60,7 @@ public: void factorizationPhaseMatrixSolver(); //Solve phase of matrix solver with a specified number of right-hand side - void solvePhaseMatrixSolver( double* solution, const int iRhsStart ,const int nRhs ); + void solvePhaseMatrixSolver( double* solution, const long long iRhsStart ,const int nRhs ); //Solve phase of matrix solver void solvePhaseMatrixSolver( double* solution ); diff --git a/src/Forward3D.cpp b/src/Forward3D.cpp index 8b8a329..0af332c 100644 --- a/src/Forward3D.cpp +++ b/src/Forward3D.cpp @@ -189,13 +189,13 @@ void Forward3D::solvePhaseForRhsConsistingInterpolatorVectors( const int numInte const int numOfEquationFinallySolved = getNumOfEquationFinallySolved(); const int numRHSDividedWithoutOdds = numInterpolatorVectors / numDivRhs; const int numAdds = numInterpolatorVectors % numDivRhs; - int iRhsStart = 0; + long long iRhsStart = 0; for( int iDiv = 0; iDiv < numDivRhs; ++iDiv ){ const int numRHSDividedWithout = iDiv < numAdds ? numRHSDividedWithoutOdds + 1 : numRHSDividedWithoutOdds; OutputFiles::m_logFile << "# Solve phase is performed simultaneously for " << numRHSDividedWithout << " right-hand sides" << ptrAnalysisControl->outputElapsedTime() << std::endl; - const int long long index = static_cast(numOfEquationFinallySolved) * static_cast(iRhsStart); + const int long long index = static_cast(numOfEquationFinallySolved) * iRhsStart; m_matrix3DAnalysis.solvePhaseMatrixSolver( &solutionForInterpolatorVectors[index], iRhsStart, numRHSDividedWithout ); - iRhsStart += numRHSDividedWithout; + iRhsStart += static_cast(numRHSDividedWithout); } //----- debug >>>>> diff --git a/src/InversionGaussNewtonDataSpace.cpp b/src/InversionGaussNewtonDataSpace.cpp index a632427..de24221 100644 --- a/src/InversionGaussNewtonDataSpace.cpp +++ b/src/InversionGaussNewtonDataSpace.cpp @@ -73,7 +73,7 @@ void InversionGaussNewtonDataSpace::inversionCalculation(){ // Perform inversion by the new method void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ - const bool useBLAS = false; + const bool useBLAS = true; // Get process ID and total process number const AnalysisControl* const ptrAnalysisControl = AnalysisControl::getInstance(); @@ -254,11 +254,11 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ const int numRHSDividedWithoutOdds = numDataThisFreq / numDivRhs; const int numAdds = numDataThisFreq % numDivRhs; - int iRhsStart = 0; + long long iRhsStart = 0; for( int iDiv = 0; iDiv < numDivRhs; ++iDiv ){ 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) * static_cast(iRhsStart); + const long long index = static_cast(numModel) * iRhsStart; transposedConstrainingMatrix.solvePhaseMatrixSolver( numRHSDivided, &sensitivityMatrix[index], &sensitivityMatrixTemp[index] );// Solve iRhsStart += numRHSDivided; } @@ -564,11 +564,11 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ CBLAS_TRANSPOSE transA = CblasNoTrans; CBLAS_TRANSPOSE transB = CblasTrans; cblas_dgemm(order, transA, transB, m, n, k, alpha, sensitivityMatrixLeft, lda, sensitivityMatrixRight, ldb, beta, result, ldc); - for( int irow = 0; irow < numDataThisFreqLeft; ++irow ){ - const int row = irow + offsetRows; - for( int icol = 0; icol < numDataThisFreqRight; ++icol ){ - const int col = icol + offsetCols; - matrix[ row * numDataTotal + col ] = result[ irow * numDataThisFreqRight + icol ]; + for( long long irow = 0; irow < numDataThisFreqLeft; ++irow ){ + const long long row = irow + static_cast(offsetRows); + for( long long icol = 0; icol < numDataThisFreqRight; ++icol ){ + const long long col = icol + static_cast(offsetCols); + matrix[ row * static_cast(numDataTotal) + col ] = result[ irow * static_cast(numDataThisFreqRight) + icol ]; } } delete [] result; @@ -897,7 +897,7 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethod() const{ // Perform inversion by the new method using inverse of [R]T[R] matrix void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMatrix() const{ - const bool useBLAS = false; + const bool useBLAS = true; // Get process ID and total process number const AnalysisControl* const ptrAnalysisControl = AnalysisControl::getInstance(); @@ -1088,11 +1088,11 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa const int numRHSDividedWithoutOdds = numDataThisFreq / numDivRhs; const int numAdds = numDataThisFreq % numDivRhs; - int iRhsStart = 0; + long long iRhsStart = 0; for( int iDiv = 0; iDiv < numDivRhs; ++iDiv ){ 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) * static_cast(iRhsStart); + const long long index = static_cast(numModel) * iRhsStart; RTRMatrix.solvePhaseMatrixSolver( numRHSDivided, &sensitivityMatrix[index], &sensitivityMatrixTemp[index] );// Solve iRhsStart += numRHSDivided; } @@ -1374,11 +1374,11 @@ void InversionGaussNewtonDataSpace::inversionCalculationByNewMethodUsingInvRTRMa CBLAS_TRANSPOSE transA = CblasNoTrans; CBLAS_TRANSPOSE transB = CblasTrans; cblas_dgemm(order, transA, transB, m, n, k, alpha, sensitivityMatrixLeft, lda, sensitivityMatrixRight, ldb, beta, result, ldc); - for( int irow = 0; irow < numDataThisFreqLeft; ++irow ){ - const int row = irow + offsetRows; - for( int icol = 0; icol < numDataThisFreqRight; ++icol ){ - const int col = icol + offsetCols; - matrix[ row * numDataTotal + col ] = result[ irow * numDataThisFreqRight + icol ]; + for( long long irow = 0; irow < numDataThisFreqLeft; ++irow ){ + const long long row = irow + static_cast(offsetRows); + for( long long icol = 0; icol < numDataThisFreqRight; ++icol ){ + const long long col = icol + static_cast(offsetCols); + matrix[ row * static_cast(numDataTotal) + col ] = result[ irow * static_cast(numDataThisFreqRight) + icol ]; } } delete [] result;