org.ejml.dense.row.decomposition.chol.CholeskyDecompositionCommon_DDRM Maven / Gradle / Ivy
/*
* Copyright (c) 2009-2017, Peter Abeles. All Rights Reserved.
*
* This file is part of Efficient Java Matrix Library (EJML).
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.ejml.dense.row.decomposition.chol;
import org.ejml.data.Complex_F64;
import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.decomposition.UtilDecompositons_DDRM;
import org.ejml.interfaces.decomposition.CholeskyDecomposition_F64;
/**
*
*
* This is an abstract class for a Cholesky decomposition. It provides the solvers, but the actual
* decomposition is provided in other classes.
*
*
* @see CholeskyDecomposition_F64
* @author Peter Abeles
*/
public abstract class CholeskyDecompositionCommon_DDRM
implements CholeskyDecomposition_F64 {
// it can decompose a matrix up to this width
protected int maxWidth=-1;
// width and height of the matrix
protected int n;
// the decomposed matrix
protected DMatrixRMaj T;
protected double[] t;
// tempoary variable used by various functions
protected double vv[];
// is it a lower triangular matrix or an upper triangular matrix
protected boolean lower;
// storage for computed determinant
protected Complex_F64 det = new Complex_F64();
/**
* Specifies if a lower or upper variant should be constructed.
*
* @param lower should a lower or upper triangular matrix be used.
*/
public CholeskyDecompositionCommon_DDRM(boolean lower) {
this.lower = lower;
}
public void setExpectedMaxSize( int numRows , int numCols ) {
if( numRows != numCols ) {
throw new IllegalArgumentException("Can only decompose square matrices");
}
this.maxWidth = numCols;
this.vv = new double[maxWidth];
}
/**
* If true the decomposition was for a lower triangular matrix.
* If false it was for an upper triangular matrix.
*
* @return True if lower, false if upper.
*/
@Override
public boolean isLower() {
return lower;
}
/**
*
* Performs Choleksy decomposition on the provided matrix.
*
*
*
* If the matrix is not positive definite then this function will return
* false since it can't complete its computations. Not all errors will be
* found. This is an efficient way to check for positive definiteness.
*
* @param mat A symmetric positive definite matrix with n ≤ widthMax.
* @return True if it was able to finish the decomposition.
*/
@Override
public boolean decompose( DMatrixRMaj mat ) {
if( mat.numRows > maxWidth ) {
setExpectedMaxSize(mat.numRows,mat.numCols);
} else if( mat.numRows != mat.numCols ) {
throw new IllegalArgumentException("Must be a square matrix.");
}
n = mat.numRows;
T = mat;
t = T.data;
if(lower) {
return decomposeLower();
} else {
return decomposeUpper();
}
}
@Override
public boolean inputModified() {
return true;
}
/**
* Performs an lower triangular decomposition.
*
* @return true if the matrix was decomposed.
*/
protected abstract boolean decomposeLower();
/**
* Performs an upper triangular decomposition.
*
* @return true if the matrix was decomposed.
*/
protected abstract boolean decomposeUpper();
@Override
public DMatrixRMaj getT(DMatrixRMaj T ) {
// write the values to T
if( lower ) {
T = UtilDecompositons_DDRM.checkZerosUT(T,n,n);
for( int i = 0; i < n; i++ ) {
for( int j = 0; j <= i; j++ ) {
T.unsafe_set(i,j,this.T.unsafe_get(i,j));
}
}
} else {
T = UtilDecompositons_DDRM.checkZerosLT(T,n,n);
for( int i = 0; i < n; i++ ) {
for( int j = i; j < n; j++ ) {
T.unsafe_set(i,j,this.T.unsafe_get(i,j));
}
}
}
return T;
}
/**
* Returns the triangular matrix from the decomposition.
*
* @return A lower or upper triangular matrix.
*/
public DMatrixRMaj getT() {
return T;
}
public double[] _getVV() {
return vv;
}
@Override
public Complex_F64 computeDeterminant() {
double prod = 1;
int total = n*n;
for( int i = 0; i < total; i += n + 1 ) {
prod *= t[i];
}
det.real = prod*prod;
det.imaginary = 0;
return det;
}
}