package com.imsl.test.example.math;

import com.imsl.math.*;

/**
 * <p>
 * Solves a small linear system with user supplied inner
 * product. </p>
 *
 * A solution to a small linear system is found. The coefficient matrix is
 * stored as a full matrix and no preconditioning is used. Typically,
 * preconditioning is required to achieve convergence in a reasonable number of
 * iterations. The user supplies a function to compute the inner product and
 * norm within the Gram-Schmidt implementation.
 *
 * @see <a href="GenMinResEx2.java">Code</a>
 * @see <a href="doc-files/GenMinResEx2.out">Output</a>
 */
public class GenMinResEx2 implements GenMinRes.Function, GenMinRes.Norm {

    // Note that if the matrix were to be read in from some
    // outside source, the code to read the matrix could reside in a constructor.
    static private double a[][] = {
        {33.0, 16.0, 72.0},
        {-24.0, -10.0, -57.0},
        {18.0, -11.0, 7.0}
    };
    static private double b[] = {129.0, -96.0, 8.5};

    /**
     * Multiplies the matrix <code>a</code> and the input vector
     * <code>p</code>. The result is returned in <code>z</code>.
     *
     * @param p a <code>double</code> array with
     * <code>p.length</code>=<code>a[0].length</code>
     * @param z a <code>double</code> array
     */
    @Override
    public void amultp(double p[], double z[]) {
        double[] result;
        result = Matrix.multiply(a, p);
        System.arraycopy(result, 0, z, 0, z.length);
    }

    /**
     * Computes the inner product of <code>x</code> and <code>y</code>.
     * 
     * The inner product \( < x, y > \) = \( \sum_{i=1}^{x.length} x[i]*y[i] \), where 
     * <code>x.length</code> \( \le \) <code>y.length</code>.
     *
     * @param x a <code>double</code> array with <code>x.length</code> \( \le \) <code>y.length</code>
     * @param y a <code>double</code> array
     * @return the inner product
     */
    @Override
    public double innerproduct(double[] x, double[] y) {
        int n = x.length;
        double tmp = 0.0;
        for (int i = 0; i < n; i++) {
            tmp += x[i] * y[i];
        }
        return tmp;
    }

    /**
     * Computes the Euclidean norm of <code>x</code>.
     *
     * The Euclidean norm of a vector, \(x \) is defined as 
     * \(\sqrt{\sum_{i=1}^{n} x_i*x_i} = \sqrt{ < x , x >} \) where \(
     * < \cdot,\cdot > \) denotes the inner product operation.
     * 
     *
     * @param x a <code>double</code> array
     * @return a <code>double</code>, the Euclidean norm of x
     */
    @Override
    public double norm(double[] x) {
        int n = x.length;
        double tmp = 0.0;
        for (int i = 0; i < n; i++) {
            tmp += x[i] * x[i];
        }
        return Math.sqrt(tmp);
    }

    public static void main(String args[]) throws Exception {
        int n = 3;

        GenMinResEx2 atp = new GenMinResEx2();

        // Construct a GenMinRes object
        GenMinRes gnmnrs = new GenMinRes(n, atp);
        gnmnrs.setVectorProducts(atp);

        // Solve Ax = b
        new PrintMatrix("x").print(gnmnrs.solve(b));
    }
}