# include # include # include # include using namespace std; # include "asa006.hpp" //****************************************************************************80 void cholesky ( double a[], int n, int nn, double u[], int *nullty, int *ifault ) //****************************************************************************80 // // Purpose: // // cholesky() computes the Cholesky factorization of an SPD matrix. // // Discussion: // // For a symmetric positive definite matrix A, the Cholesky factor U // is an upper triangular matrix such that A = U' * U. // // This routine was originally named "CHOL", but that conflicted with // a built in MATLAB routine name. // // The missing initialization "II = 0" has been added to the code. // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 12 February 2008 // // Author: // // Original FORTRAN77 version by Michael Healy. // Modifications by AJ Miller. // C++ version by John Burkardt. // // Reference: // // Michael Healy, // Algorithm AS 6: // Triangular decomposition of a symmetric matrix, // Applied Statistics, // Volume 17, Number 2, 1968, pages 195-197. // // Parameters: // // Input, double A((N*(N+1))/2), a symmetric positive definite matrix // stored by rows in lower triangular form as a one dimensional array, // in the sequence // A(1,1), // A(2,1), A(2,2), // A(3,1), A(3,2), A(3,3), and so on. // // Input, int N, the order of A. // // Input, int NN, the dimension of the array used to store A, // which should be at least (N*(N+1))/2. // // Output, double U((N*(N+1))/2), an upper triangular matrix, // stored by columns, which is the Cholesky factor of A. The program is // written in such a way that A and U can share storage. // // Output, int NULLTY, the rank deficiency of A. If NULLTY is zero, // the matrix is judged to have full rank. // // Output, int IFAULT, an error indicator. // 0, no error was detected; // 1, if N < 1; // 2, if A is not positive semi-definite. // 3, if NN < (N*(N+1))/2. // // Local: // // double ETA, should be set equal to the smallest positive // value such that 1.0 + ETA is calculated as being greater than 1.0 in the // accuracy being used. // { double eta = 1.0E-09; int i; int icol; int ii; int irow; int j; int k; int kk; int l; int m; double w; double x; *ifault = 0; *nullty = 0; if ( n <= 0 ) { *ifault = 1; return; } if ( nn < ( n * ( n + 1 ) ) / 2 ) { *ifault = 3; return; } j = 1; k = 0; ii = 0; // // Factorize column by column, ICOL = column number. // for ( icol = 1; icol <= n; icol++ ) { ii = ii + icol; x = eta * eta * a[ii-1]; l = 0; kk = 0; // // IROW = row number within column ICOL. // for ( irow = 1; irow <= icol; irow++ ) { kk = kk + irow; k = k + 1; w = a[k-1]; m = j; for ( i = 1; i < irow; i++ ) { l = l + 1; w = w - u[l-1] * u[m-1]; m = m + 1; } l = l + 1; if ( irow == icol ) { break; } if ( u[l-1] != 0.0 ) { u[k-1] = w / u[l-1]; } else { u[k-1] = 0.0; if ( fabs ( x * a[k-1] ) < w * w ) { *ifault = 2; return; } } } // // End of row, estimate relative accuracy of diagonal element. // if ( fabs ( w ) <= fabs ( eta * a[k-1] ) ) { u[k-1] = 0.0; *nullty = *nullty + 1; } else { if ( w < 0.0 ) { *ifault = 2; return; } u[k-1] = sqrt ( w ); } j = j + icol; } return; } //****************************************************************************80 void subchl ( double a[], int b[], int n, double u[], int *nullty, int *ifault, int ndim, double *det ) //****************************************************************************80 // // Purpose: // // subchl() computes the Cholesky factorization of a subset of an SPD matrix. // // Licensing: // // This code is distributed under the MIT license. // // Modified: // // 11 February 2008 // // Author: // // Original FORTRAN77 version by Michael Healy, PR Freeman. // C++ version by John Burkardt. // // Reference: // // PR Freeman, // Remark AS R44: // A Remark on AS 6 and AS7: Triangular decomposition of a symmetric matrix // and Inversion of a positive semi-definite symmetric matrix, // Applied Statistics, // Volume 31, Number 3, 1982, pages 336-339. // // Michael Healy, // Algorithm AS 6: // Triangular decomposition of a symmetric matrix, // Applied Statistics, // Volume 17, Number 2, 1968, pages 195-197. // // Parameters: // // Input, double A((M*(M+1))/2), a symmetric positive definite matrix // stored by rows in lower triangular form as a one dimensional array, // in the sequence // A(1,1), // A(2,1), A(2,2), // A(3,1), A(3,2), A(3,3), and so on. // In the simplest case, M, the order of A, is equal to N. // // Input, int B(N), indicates the order in which the // rows and columns of A are to be used. In the simplest case, // B = (1,2,3...,N). // // Input, int N, the order of the matrix, that is, // the matrix formed by using B to select N rows and columns of A. // // Output, double U((N*(N+1))/2), an upper triangular matrix, // stored by columns, which is the Cholesky factor of A. The program is // written in such a way that A and U can share storage. // // Output, int *NULLTY, the rank deficiency of A. // If NULLTY is zero, the matrix is judged to have full rank. // // Output, int *IFAULT, an error indicator. // 0, no error was detected; // 1, if N < 1; // 2, if A is not positive semi-definite. // // Input, int NDIM, the dimension of A and U, which might // be presumed to be (N*(N+1))/2. // // Output, double *DET, the determinant of the matrix. // { double eta = 1.0E-09; int i; int icol; int ii; int ij; int irow; int j; int jj; int k; int kk; int l; int m; double w; double x; *ifault = 0; *nullty = 0; *det = 1.0; if ( n <= 0 ) { *ifault = 1; return; } *ifault = 2; j = 1; k = 0; for ( icol = 1; icol <= n; icol++ ) { ij = ( b[icol-1] * ( b[icol-1] - 1 ) ) / 2; ii = ij + b[icol-1]; x = eta * eta * a[ii-1]; l = 0; for ( irow = 1; irow <= icol; irow++ ) { kk = ( b[irow-1] * ( b[irow-1] + 1 ) ) / 2; k = k + 1; jj = ij + b[irow-1]; w = a[jj-1]; m = j; for ( i = 1; i <= irow - 1; i++ ) { l = l + 1; w = w - u[l-1] * u[m-1]; m = m + 1; } l = l + 1; if ( irow == icol ) { break; } if ( u[l-1] != 0.0 ) { u[k-1] = w / u[l-1]; } else { if ( fabs ( x * a[kk-1] ) < w * w ) { *ifault = 2; return; } u[k-1] = 0.0; } } if ( fabs ( eta * a[kk-1] ) <= fabs ( w ) ) { if ( w < 0.0 ) { *ifault = 2; return; } u[k-1] = sqrt ( w ); } else { u[k-1] = 0.0; *nullty = *nullty + 1; } j = j + icol; *det = *det * u[k-1] * u[k-1]; } return; }