001    /*
002     * Copyright (C) 2008-2010 by Holger Arndt
003     *
004     * This file is part of the Universal Java Matrix Package (UJMP).
005     * See the NOTICE file distributed with this work for additional
006     * information regarding copyright ownership and licensing.
007     *
008     * UJMP is free software; you can redistribute it and/or modify
009     * it under the terms of the GNU Lesser General Public License as
010     * published by the Free Software Foundation; either version 2
011     * of the License, or (at your option) any later version.
012     *
013     * UJMP is distributed in the hope that it will be useful,
014     * but WITHOUT ANY WARRANTY; without even the implied warranty of
015     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
016     * GNU Lesser General Public License for more details.
017     *
018     * You should have received a copy of the GNU Lesser General Public
019     * License along with UJMP; if not, write to the
020     * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
021     * Boston, MA  02110-1301  USA
022     */
023    
024    package org.ujmp.core.doublematrix.calculation.general.decomposition;
025    
026    import org.ujmp.core.Matrix;
027    import org.ujmp.core.doublematrix.DenseDoubleMatrix2D;
028    import org.ujmp.core.exceptions.MatrixException;
029    import org.ujmp.core.util.DecompositionOps;
030    import org.ujmp.core.util.UJMPSettings;
031    
032    public interface InvSymm<T> extends Inv<T> {
033    
034            public static int THRESHOLD = 100;
035    
036            public T calc(T source);
037    
038            public static final InvSymm<Matrix> MATRIX = new InvSymm<Matrix>() {
039    
040                    public final Matrix calc(Matrix source) {
041                            if (source.getDimensionCount() != 2 || source.getRowCount() != source.getColumnCount()) {
042                                    throw new MatrixException(
043                                                    "inverse only possible for square matrices. use pinv or ginv instead");
044                            }
045                            if (UJMPSettings.getNumberOfThreads() == 1) {
046                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
047                                            return MATRIXLARGESINGLETHREADED.calc(source);
048                                    } else {
049                                            return MATRIXSMALLSINGLETHREADED.calc(source);
050                                    }
051                            } else {
052                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
053                                            return MATRIXLARGEMULTITHREADED.calc(source);
054                                    } else {
055                                            return MATRIXSMALLMULTITHREADED.calc(source);
056                                    }
057                            }
058                    }
059            };
060    
061            public final InvSymm<Matrix> UJMP = new InvSymm<Matrix>() {
062                    public final Matrix calc(Matrix source) {
063                            DenseDoubleMatrix2D b = DenseDoubleMatrix2D.factory.zeros(source.getRowCount(), source
064                                            .getRowCount());
065                            for (int i = (int) source.getRowCount(); --i >= 0;) {
066                                    b.setDouble(1.0, i, i);
067                            }
068                            return LU.INSTANCE.solve(source, b);
069                    }
070            };
071    
072            public final InvSymm<Matrix> INSTANCE = MATRIX;
073    
074            public static final InvSymm<Matrix> MATRIXSMALLSINGLETHREADED = UJMP;
075    
076            public static final InvSymm<Matrix> MATRIXLARGESINGLETHREADED = new InvSymm<Matrix>() {
077                    public final Matrix calc(Matrix source) {
078                            // no special implementation for symmetric matrices
079                            Inv<Matrix> inv = null;
080                            if (UJMPSettings.isUseJBlas()) {
081                                    inv = DecompositionOps.INV_JBLAS;
082                            }
083                            if (inv == null && UJMPSettings.isUseOjalgo()) {
084                                    inv = DecompositionOps.INV_OJALGO;
085                            }
086                            if (inv == null && UJMPSettings.isUseEJML()) {
087                                    inv = DecompositionOps.INV_EJML;
088                            }
089                            if (inv == null && UJMPSettings.isUseMTJ()) {
090                                    inv = DecompositionOps.INV_MTJ;
091                            }
092                            if (inv == null) {
093                                    inv = UJMP;
094                            }
095                            return inv.calc(source);
096                    }
097            };
098    
099            public static final InvSymm<Matrix> MATRIXLARGEMULTITHREADED = new InvSymm<Matrix>() {
100                    public Matrix calc(Matrix source) {
101                            // no special implementation for symmetric matrices
102                            Inv<Matrix> inv = null;
103                            if (UJMPSettings.isUseJBlas()) {
104                                    inv = DecompositionOps.INV_JBLAS;
105                            }
106                            if (inv == null && UJMPSettings.isUseOjalgo()) {
107                                    inv = DecompositionOps.INV_OJALGO;
108                            }
109                            if (inv == null && UJMPSettings.isUseEJML()) {
110                                    inv = DecompositionOps.INV_EJML;
111                            }
112                            if (inv == null && UJMPSettings.isUseMTJ()) {
113                                    inv = DecompositionOps.INV_MTJ;
114                            }
115                            if (inv == null) {
116                                    inv = UJMP;
117                            }
118                            return inv.calc(source);
119                    }
120            };
121    
122            public static final InvSymm<Matrix> MATRIXSMALLMULTITHREADED = UJMP;
123    
124    }