/**
 * Copyright (C) 2011 - present by OpenGamma Inc. and the OpenGamma group of companies
 * 
 * Please see distribution for license.
 */
package com.opengamma.analytics.math.linearalgebra;

import org.apache.commons.lang.Validate;

import com.opengamma.analytics.math.matrix.DoubleMatrix1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix2D;
import com.opengamma.analytics.math.matrix.MatrixAlgebra;
import com.opengamma.analytics.math.matrix.OGMatrixAlgebra;

/**
 * Results of the OpenGamma implementation of Cholesky decomposition.
 */
public class CholeskyDecompositionOpenGammaResult implements CholeskyDecompositionResult {

  private static final MatrixAlgebra ALGEBRA = new OGMatrixAlgebra();
  /**
   * The array that store the data.
   */
  private final double[][] _lArray;
  /**
   * The matrix L, result of the decomposition.
   */
  private final DoubleMatrix2D _l;
  /**
   * The matrix L^T, result of the decomposition.
   */
  private final DoubleMatrix2D _lT;
  /**
   * The determinant of the original matrix A = L L^T.
   */
  private double _determinant;

  /**
   * Constructor.
   * @param lArray The matrix L as an array of doubles.
   */
  public CholeskyDecompositionOpenGammaResult(final double[][] lArray) {
    _lArray = lArray;
    _l = new DoubleMatrix2D(_lArray);
    _lT = ALGEBRA.getTranspose(_l);
    _determinant = 1.0;
    for (int loopdiag = 0; loopdiag < _lArray.length; ++loopdiag) {
      _determinant *= _lArray[loopdiag][loopdiag] * _lArray[loopdiag][loopdiag];
    }
  }

  @Override
  public DoubleMatrix1D solve(DoubleMatrix1D b) {
    return new DoubleMatrix1D(b.getData());
  }

  @Override
  public double[] solve(double[] b) {
    int dim = b.length;
    Validate.isTrue(dim == _lArray.length, "b array of incorrect size");
    final double[] x = new double[dim];
    System.arraycopy(b, 0, x, 0, dim);
    // L y = b (y stored in x array)
    for (int looprow = 0; looprow < dim; looprow++) {
      x[looprow] /= _lArray[looprow][looprow];
      for (int j = looprow + 1; j < dim; j++) {
        x[j] -= x[looprow] * _lArray[j][looprow];
      }
    }
    // L^T x = y
    for (int looprow = dim - 1; looprow >= -0; looprow--) {
      x[looprow] /= _lArray[looprow][looprow];
      for (int j = 0; j < looprow; j++) {
        x[j] -= x[looprow] * _lArray[looprow][j];
      }
    }
    return x;
  }

  @Override
  public DoubleMatrix2D solve(DoubleMatrix2D b) {
    int nbRow = b.getNumberOfRows();
    int nbCol = b.getNumberOfColumns();
    Validate.isTrue(nbRow == _lArray.length, "b array of incorrect size");
    double[][] bArray = b.getData();
    final double[][] x = new double[nbRow][nbCol];
    for (int looprow = 0; looprow < nbRow; looprow++) {
      System.arraycopy(bArray[looprow], 0, x[looprow], 0, nbCol);
    }
    // L Y = B (Y stored in x array)
    for (int loopcol = 0; loopcol < nbCol; loopcol++) {
      for (int looprow = 0; looprow < nbRow; looprow++) {
        x[looprow][loopcol] /= _lArray[looprow][looprow];
        for (int j = looprow + 1; j < nbRow; j++) {
          x[j][loopcol] -= x[looprow][loopcol] * _lArray[j][looprow];
        }
      }
    }
    // L^T X = Y
    for (int loopcol = 0; loopcol < nbCol; loopcol++) {
      for (int looprow = nbRow - 1; looprow >= -0; looprow--) {
        x[looprow][loopcol] /= _lArray[looprow][looprow];
        for (int j = 0; j < looprow; j++) {
          x[j][loopcol] -= x[looprow][loopcol] * _lArray[looprow][j];
        }
      }
    }
    return new DoubleMatrix2D(x);
  }

  @Override
  public DoubleMatrix2D getL() {
    return _l;
  }

  @Override
  public DoubleMatrix2D getLT() {
    return _lT;
  }

  @Override
  public double getDeterminant() {
    return _determinant;
  }

}
