/****************************************************************************** * * AMDiS - Adaptive multidimensional simulations * * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved. * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis * * Authors: * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al. * * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. * * * This file is part of AMDiS * * See also license.opensource.txt in the distribution. * ******************************************************************************/ #include "Cholesky.h" bool Cholesky::factorization(Matrix *A, Vector *p) { FUNCNAME("Cholesky::factorization()"); int n = A->getNumRows(); // Checking memory for vector P of diagonal elements of factorization. static Vector *pT = nullptr; if (p) { if (p->getSize() != n) p->resize(n); } else { if (pT) { if (pT->getSize() != n) pT->resize(n); } else pT = new Vector(n); p = pT; } // Constructs the Cholesky factorization A = L*L, with L lower triangular. // Only the upper triangle need be given; it is not modified. // The Cholesky factor L is stored in the lower triangle of A, except for // its diagonal which is stored in P. int i, j, k; double sum; for (i=0; i=0; k--) sum -= (*A)[i][k] * (*A)[j][k]; if (i==j) { if (sum<=0) { MSG("Matrix is (numerically) not positive definite!\n"); MSG("Cholesky decomposition does not work; choose another method for solving the system.\n"); return false; } (*p)[i] = sqrt(sum); } else (*A)[j][i] = sum / (*p)[i]; } } return true; } bool Cholesky::solve(Matrix *A, Vector *b, Vector *x, Vector *p) { FUNCNAME("Cholesky::solve"); bool success = true; int n = b->getSize(); TEST_EXIT(n == A->getNumRows()) ("Dimensions of matrix and vector do not match!\n"); // Checking memory for solution vector X. if (x && (x->getSize() != n)) x->resize(n); if (!x) x = new Vector(n); // Checking vector P. static Vector *pT = nullptr; if (!p || (p->getSize() != n)) { if (pT && pT->getSize() != n) delete pT; if (!pT) pT = new Vector(n); p = pT; success = factorization(A, p); } // Now solve the system by backsubstitution. int i, k; double sum; for (i=0; i=0; k--) sum -= (*A)[i][k] * (*x)[k]; (*x)[i] = sum / (*p)[i]; } for (i=n-1; i>=0; i--) // Solve L^T*X = Y. { for (sum=(*x)[i], k=i+1; k *A, Vector > *b, Vector > *x, Vector *p) { FUNCNAME("Cholesky::solve"); bool success = true; int n = b->getSize(); TEST_EXIT(n == A->getNumRows()) ("Dimensions of matrix and vector do not match!\n"); // Checking memory for solution vector X. if (x && (x->getSize() != n)) x->resize(n); if (!x) x = new Vector >(n); // Checking vector P. static Vector *pT = nullptr; if (!p || (p->getSize() != n)) { pT = new Vector(n); p = pT; success = factorization(A, p); } // Now solve the system by backsubstitution. int i, k; WorldVector vec_sum; for (i=0; i=0; k--) vec_sum -= (*x)[k] * (*A)[i][k] ; (*x)[i] = vec_sum * (1.0/(*p)[i]); } for (i=n-1; i>=0; i--) // Solve L^T*X = Y. { for (vec_sum=(*x)[i], k=i+1; k