org.ejml.dense.block.linsol.chol.CholeskyOuterSolver_DDRB 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.block.linsol.chol;

import org.ejml.data.DMatrixRBlock;
import org.ejml.data.DSubmatrixD1;
import org.ejml.dense.block.MatrixOps_DDRB;
import org.ejml.dense.block.TriangularSolver_DDRB;
import org.ejml.dense.block.decomposition.chol.CholeskyOuterForm_DDRB;
import org.ejml.dense.row.SpecializedOps_DDRM;
import org.ejml.interfaces.decomposition.CholeskyDecomposition_F64;
import org.ejml.interfaces.linsol.LinearSolverDense;


/**
 *  Linear solver that uses a block cholesky decomposition. 
 *
 * 
 * Solver works by using the standard Cholesky solving strategy:

 * A=L*L^T 

 * A*x=b

 * L*L^T*x = b 

 * L*y = b

 * L^T*x = y

 * x = L^-Ty
 * 
 *
 * 
 * It is also possible to use the upper triangular cholesky decomposition.
 * 
 *
 * @author Peter Abeles
 */
public class CholeskyOuterSolver_DDRB implements LinearSolverDense {

    // cholesky decomposition
    private CholeskyOuterForm_DDRB decomposer = new CholeskyOuterForm_DDRB(true);

    // size of a block take from input matrix
    private int blockLength;

    // temporary data structure used in some calculation.
    private double temp[];

    /**
     * Decomposes and overwrites the input matrix.
     *
     * @param A Semi-Positive Definite (SPD) system matrix. Modified. Reference saved.
     * @return If the matrix can be decomposed.  Will always return false of not SPD.
     */
    @Override
    public boolean setA(DMatrixRBlock A) {
        // Extract a lower triangular solution
        if( !decomposer.decompose(A) )
            return false;

        blockLength = A.blockLength;

        return true;
    }

    @Override
    public /**/double quality() {
        return SpecializedOps_DDRM.qualityTriangular(decomposer.getT(null));
    }

    /**
     * If X == null then the solution is written into B.  Otherwise the solution is copied
     * from B into X.
     */
    @Override
    public void solve(DMatrixRBlock B, DMatrixRBlock X) {
        if( B.blockLength != blockLength )
            throw new IllegalArgumentException("Unexpected blocklength in B.");

        DSubmatrixD1 L = new DSubmatrixD1(decomposer.getT(null));

        if( X != null ) {
            if( X.blockLength != blockLength )
                throw new IllegalArgumentException("Unexpected blocklength in X.");
            if( X.numRows != L.col1 ) throw new IllegalArgumentException("Not enough rows in X");
        }
        
        if( B.numRows != L.col1 ) throw new IllegalArgumentException("Not enough rows in B");

        //  L * L^T*X = B

        // Solve for Y:  L*Y = B
        TriangularSolver_DDRB.solve(blockLength,false,L,new DSubmatrixD1(B),false);

        // L^T * X = Y
        TriangularSolver_DDRB.solve(blockLength,false,L,new DSubmatrixD1(B),true);

        if( X != null ) {
            // copy the solution from B into X
            MatrixOps_DDRB.extractAligned(B,X);
        }

    }

    @Override
    public void invert(DMatrixRBlock A_inv) {
        DMatrixRBlock T = decomposer.getT(null);
        if( A_inv.numRows != T.numRows || A_inv.numCols != T.numCols )
            throw new IllegalArgumentException("Unexpected number or rows and/or columns");


        if( temp == null || temp.length < blockLength*blockLength )
            temp = new double[ blockLength* blockLength ];

        // zero the upper triangular portion of A_inv
        MatrixOps_DDRB.zeroTriangle(true,A_inv);

        DSubmatrixD1 L = new DSubmatrixD1(T);
        DSubmatrixD1 B = new DSubmatrixD1(A_inv);

        // invert L from cholesky decomposition and write the solution into the lower
        // triangular portion of A_inv
        // B = inv(L)
        TriangularSolver_DDRB.invert(blockLength,false,L,B,temp);

        // B = L^-T * B
        // todo could speed up by taking advantage of B being lower triangular
        // todo take advantage of symmetry
        TriangularSolver_DDRB.solveL(blockLength,L,B,true);
    }

    @Override
    public boolean modifiesA() {
        return decomposer.inputModified();
    }

    @Override
    public boolean modifiesB() {
        return true;
    }

    @Override
    public CholeskyDecomposition_F64 getDecomposition() {
        return decomposer;
    }
}