From 266d9f87baa52498b6dcae10b7cf5d6afd0ad27b Mon Sep 17 00:00:00 2001 From: qichenghua Date: Wed, 19 Apr 2017 10:52:39 +0200 Subject: [PATCH 1/2] added new linear solver based on singular value decomposition in Nox --- .../nox/src-lapack/NOX_LAPACK_LinearSolver.H | 34 ++++++++++++++++--- 1 file changed, 29 insertions(+), 5 deletions(-) diff --git a/trilinos-nox/packages/nox/src-lapack/NOX_LAPACK_LinearSolver.H b/trilinos-nox/packages/nox/src-lapack/NOX_LAPACK_LinearSolver.H index 6a009e83d..2b031c8b6 100644 --- a/trilinos-nox/packages/nox/src-lapack/NOX_LAPACK_LinearSolver.H +++ b/trilinos-nox/packages/nox/src-lapack/NOX_LAPACK_LinearSolver.H @@ -136,6 +136,13 @@ namespace NOX { //! LAPACK wrappers Teuchos::LAPACK lapack; + + + double* sv;//singular values + int rcond; + int rank; + double* work; + int lwork; }; @@ -150,7 +157,9 @@ NOX::LAPACK::LinearSolver::LinearSolver(int n) : pivots(n), isValidLU(false), blas(), - lapack() + lapack(), + sv(NULL), + work(NULL) { } @@ -161,13 +170,17 @@ NOX::LAPACK::LinearSolver::LinearSolver(const NOX::LAPACK::LinearSolver& s pivots(s.pivots), isValidLU(s.isValidLU), blas(), - lapack() + lapack(), + sv(NULL), + work(NULL) { } template NOX::LAPACK::LinearSolver::~LinearSolver() { + if(sv) delete [] sv; + if(work) delete [] work; } template @@ -179,6 +192,8 @@ NOX::LAPACK::LinearSolver::operator=(const NOX::LAPACK::LinearSolver& s) lu = s.lu; pivots = s.pivots; isValidLU = s.isValidLU; + sv = NULL;//sv=((s.sv) ? (new double[mat.numRows()]) : NULL); + work = NULL; } return *this; @@ -234,8 +249,17 @@ NOX::LAPACK::LinearSolver::solve(bool trans, int ncols, T* output) if (!isValidLU) { lu = mat; lapack.GETRF(n, n, &lu(0,0), n, &pivots[0], &info); - if (info != 0) - return false; + if (info != 0){//try solving a least squares problem instead (in case the linear system has infinitely many solutions) + rcond = -1.0; + lwork = 10*n+1;//max(5*n,1) is the minimum, the more, the better (But the more, the less memory efficient) + if (sv) delete [] sv; + if (work) delete [] work; + sv=new double [n]; + work=new double [lwork]; + lu=mat; + lapack.GELSS(n,n,ncols,&lu(0,0), n, output, n, sv, rcond, &rank, work, lwork, &info); + return (info == 0); + } isValidLU = true; } @@ -250,7 +274,7 @@ NOX::LAPACK::LinearSolver::solve(bool trans, int ncols, T* output) lapack.GETRS(tr, n, ncols, &lu(0,0), n, &pivots[0], output, n, &info); if (info != 0) - return false; + return false;//this throws an error in the higher level function NOX::Direction::Newton::compute or prints a Warning there. return true; } From a2358fc90e9a19be55f1adba2d5592c82e6b35fa Mon Sep 17 00:00:00 2001 From: qichenghua Date: Wed, 26 Apr 2017 18:06:54 +0200 Subject: [PATCH 2/2] fixed minor bugs in NOX_LAPACK_LINEARSOLVER_H --- .../nox/src-lapack/NOX_LAPACK_LinearSolver.H | 31 +++++++++---------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/trilinos-nox/packages/nox/src-lapack/NOX_LAPACK_LinearSolver.H b/trilinos-nox/packages/nox/src-lapack/NOX_LAPACK_LinearSolver.H index 2b031c8b6..f01af1903 100644 --- a/trilinos-nox/packages/nox/src-lapack/NOX_LAPACK_LinearSolver.H +++ b/trilinos-nox/packages/nox/src-lapack/NOX_LAPACK_LinearSolver.H @@ -138,11 +138,10 @@ namespace NOX { Teuchos::LAPACK lapack; - double* sv;//singular values - int rcond; + int lwork; int rank; + double* sv;//singular values double* work; - int lwork; }; @@ -158,8 +157,9 @@ NOX::LAPACK::LinearSolver::LinearSolver(int n) : isValidLU(false), blas(), lapack(), - sv(NULL), - work(NULL) + lwork(10*mat.numRows()+1), + sv(new double [mat.numRows()]),//this array is only used as output and is not necessary + work(new double [lwork])//this array is only used as output and is not necessary { } @@ -171,8 +171,9 @@ NOX::LAPACK::LinearSolver::LinearSolver(const NOX::LAPACK::LinearSolver& s isValidLU(s.isValidLU), blas(), lapack(), - sv(NULL), - work(NULL) + lwork(10*mat.numRows()+1), + sv(new double [mat.numRows()]),//this array is only used as output and is not necessary + work(new double [lwork])//this array is only used as output and is not necessary { } @@ -192,8 +193,9 @@ NOX::LAPACK::LinearSolver::operator=(const NOX::LAPACK::LinearSolver& s) lu = s.lu; pivots = s.pivots; isValidLU = s.isValidLU; - sv = NULL;//sv=((s.sv) ? (new double[mat.numRows()]) : NULL); - work = NULL; + lwork = s.lwork; + sv = new double [mat.numRows()];//this array is only used as output and is not necessary + work = new double [lwork];//this array is only used as output and is not necessary } return *this; @@ -244,21 +246,16 @@ NOX::LAPACK::LinearSolver::solve(bool trans, int ncols, T* output) { int info; int n = mat.numRows(); + int rcond = -1.0; // Compute LU factorization if necessary if (!isValidLU) { lu = mat; lapack.GETRF(n, n, &lu(0,0), n, &pivots[0], &info); if (info != 0){//try solving a least squares problem instead (in case the linear system has infinitely many solutions) - rcond = -1.0; - lwork = 10*n+1;//max(5*n,1) is the minimum, the more, the better (But the more, the less memory efficient) - if (sv) delete [] sv; - if (work) delete [] work; - sv=new double [n]; - work=new double [lwork]; lu=mat; lapack.GELSS(n,n,ncols,&lu(0,0), n, output, n, sv, rcond, &rank, work, lwork, &info); - return (info == 0); + return false; } isValidLU = true; } @@ -274,7 +271,7 @@ NOX::LAPACK::LinearSolver::solve(bool trans, int ncols, T* output) lapack.GETRS(tr, n, ncols, &lu(0,0), n, &pivots[0], output, n, &info); if (info != 0) - return false;//this throws an error in the higher level function NOX::Direction::Newton::compute or prints a Warning there. + return false; return true; }