001/*
002 *                    BioJava development code
003 *
004 * This code may be freely distributed and modified under the
005 * terms of the GNU Lesser General Public Licence.  This should
006 * be distributed with the code.  If you do not have a copy,
007 * see:
008 *
009 *      http://www.gnu.org/copyleft/lesser.html
010 *
011 * Copyright for this code is held jointly by the individual
012 * authors.  These should be listed in @author doc comments.
013 *
014 * For more information on the BioJava project and its aims,
015 * or to join the biojava-l mailing list, visit the home page
016 * at:
017 *
018 *      http://www.biojava.org/
019 *
020 */
021package org.biojava.nbio.survival.cox.stats;
022
023/**
024 *
025 * @author Scooter Willis <willishf at gmail dot com>
026 */
027public class Cholesky2 {
028
029           /* $Id: cholesky2.c 11357 2009-09-04 15:22:46Z therneau $
030         **
031         ** subroutine to do Cholesky decompostion on a matrix: C = FDF'
032         **   where F is lower triangular with 1's on the diagonal, and D is diagonal
033         **
034         ** arguments are:
035         **     n         the size of the matrix to be factored
036         **     **matrix  a ragged array containing an n by n submatrix to be factored
037         **     toler     the threshold value for detecting "singularity"
038         **
039         **  The factorization is returned in the lower triangle, D occupies the
040         **    diagonal and the upper triangle is left undisturbed.
041         **    The lower triangle need not be filled in at the start.
042         **
043         **  Return value:  the rank of the matrix (non-negative definite), or -rank
044         **     it not SPD or NND
045         **
046         **  If a column is deemed to be redundant, then that diagonal is set to zero.
047         **
048         **   Terry Therneau
049         */
050        /**
051         *
052         * @param matrix
053         * @param n
054         * @param toler
055         * @return
056         */
057        public static int process(double[][] matrix, int n, double toler) {
058                double temp;
059                int i, j, k;
060                double eps, pivot;
061                int rank;
062                int nonneg;
063
064                nonneg = 1;
065                eps = 0;
066                for (i = 0; i < n; i++) {
067                        if (matrix[i][i] > eps) {
068                                eps = matrix[i][i];
069                        }
070                        for (j = (i + 1); j < n; j++) {
071                                matrix[j][i] = matrix[i][j];
072                        }
073                }
074                eps *= toler;
075
076                rank = 0;
077                for (i = 0; i < n; i++) {
078                        pivot = matrix[i][i];
079                        if (pivot < eps) {
080                                matrix[i][i] = 0;
081                                if (pivot < -8 * eps) {
082                                        nonneg = -1;
083                                }
084                        } else {
085                                rank++;
086                                for (j = (i + 1); j < n; j++) {
087                                        temp = matrix[j][i] / pivot;
088                                        matrix[j][i] = temp;
089                                        matrix[j][j] -= temp * temp * pivot;
090                                        for (k = (j + 1); k < n; k++) {
091                                                matrix[k][j] -= temp * matrix[k][i];
092                                        }
093                                }
094                        }
095                }
096                return (rank * nonneg);
097        }
098
099}