/****************************************************************************
 *
 * Bidiagonal.c
 *
 * Abstract
 *   To implement the bidiagonalization of general complex matrices by
 *   invoking LAPACK's routine. Bidiagonalization of a complex matrix
 *   is used to compare the performance with the block Lanczos
 *   tridiagonalization algorithm.
 *
 * Student:  Guohong Liu
 *
 * ID:       0385117
 *
 ****************************************************************************/

#include <stdio.h>
#include <malloc.h>
#include "bidiagonal.h"
#include "blklan.h"

/***************************************************************************
 *
 * Abstract:
 *   Reducing a general complex matrix A to bidiagonal form:
 *      A = Q * B * P'.
 *
 * Input
 *   n      Order of matrix A.
 *   A_m    Pointer to the memory of matrix A.
 *
 * Output
 *   d_m    Pointer to the memory of vector d.
 *   b_m    Pointer to the memory of vector b.
 *   Q_m    Pointer to the memory of matrix Q.
 *   P_m    Pointer to the memory of matrix P.
 *
 * Return
 *   0      Successful exit
 *   1      Memory allocation exception occurs
 *
 * NOTE
 *   P in fact is output as P', not P
 *
 ***************************************************************************/

int brd (int n, DoubleComplex *A_m, DoubleReal *d_m, DoubleReal *b_m,
         DoubleComplex *Q_m, DoubleComplex *P_m)
{
    DoubleComplexMat A, Q, P, R, A1;
    DoubleRealVec d, b;
    DoubleComplex *tauq, *taup, *work;
    int info;

    A.m = n;
    A.n = n;
    A.mat = A_m;
    Q.m = n;
    Q.n = n;
    Q.mat = Q_m;
    P.m = n;
    P.n = n;
    P.mat = P_m;
    d.n = n;
    d.vec = d_m;
    b.n = n - 1;
    b.vec = b_m;

    /* Matrix A will be overwriten in bidiagonal reduction,
     *  so copy it to A1.
     */
    A1.m = n;
    A1.n = n;
    if (!(A1.mat = (DoubleComplex *)
                    malloc (n * n * sizeof (DoubleComplex)))) {
        printf ("brd: Malloc for A1.mat failed.\n");
        return 1;
    }
    zlacpy_ ("", &n, &n, A.mat, &n, A1.mat, &n);


    /* Allocate memory for *tauq, *taup, *work   */
    if (!(tauq = (DoubleComplex *)malloc (n * sizeof (DoubleComplex)))) {
        printf ("brd: Malloc for *tauq failed.\n");
        return 1;
    }
    if (!(taup = (DoubleComplex *)malloc (n * sizeof (DoubleComplex)))) {
        printf ("brd: Malloc for *taup failed.\n");
        return 1;
    }
    if (!(work = (DoubleComplex *)malloc (n * sizeof (DoubleComplex)))) {
        printf ("brd: Malloc for *work failed.\n");
        return 1;
    }

    zgebrd_ (&n, &n, A1.mat, &n, d.vec, b.vec, tauq, taup, work, &n, &info);

    /* Copy A1.mat to Q->mat because zungbr_ will overwrite A1.mat as
     * the result of Q, or P.
     */
    zlacpy_ ("", &n, &n, A1.mat, &n, Q.mat, &n);
    zlacpy_ ("", &n, &n, A1.mat, &n, P.mat, &n);

    /* Compute Q.       */
    zungbr_ ("Q", &n, &n, &n, Q.mat, &n, tauq, work, &n, &info);
    /* It computes P', not just P. It is unlike the computation of Q   */
    zungbr_ ("P", &n, &n, &n, P.mat, &n, taup, work, &n, &info);

    /* Release memory.  */
    free (A1.mat);
    free (tauq);
    free (taup);
    free (work);
    return 0;
}

/***************************************************************************
 *
 * Abstract
 *   Check the errors in factorization and orthogonality in
 *   bidiagonalization algorith.
 *
 * Input
 *   n      Order of matrix A.
 *   A_m    Pointer to the memory of matrix A.
 *   d_m    Pointer to the memory of vector d.
 *   b_m    Pointer to the memory of vector b.
 *   Q_m    Pointer to the memory of matrix Q.
 *   P_m    Pointer to the memory of matrix P.
 *
 *
 * Output
 *   None
 *
 * Return
 *   0      Successful exit
 *   1      Memory allocation exception occurs
 *
 * NOTE
 *   P is in fact specified as P'.
 *
 ***************************************************************************/

int brdErrChk (int n, DoubleComplex *A_m, DoubleReal *d_m, DoubleReal *b_m,
               DoubleComplex *Q_m, DoubleComplex *P_m)
{
    double err;
    DoubleComplexMat A, Q, P, R, Tmp;
    DoubleRealVec d, b;
    DoubleComplex *mat;
    int i, j;

    A.m = n;
    A.n = n;
    A.mat = A_m;
    Q.m = n;
    Q.n = n;
    Q.mat = Q_m;
    P.m = n;
    P.n = n;
    P.mat = P_m;
    d.n = n;
    d.vec = d_m;
    b.n = n - 1;
    b.vec = b_m;

    orthErr (&Q, &err);
    printf ("\nError in orthogonality of Q: %1.3e\n", err);
    orthErr (&P, &err);
    printf ("Error in orthogonality of P: %1.3e\n", err);

    /* Q' * A * (P')' - B   */
    if (!(R.mat = (DoubleComplex *)
                    malloc (n * n * sizeof (DoubleComplex)))) {
        printf ("brdFactErr: Malloc for R.mat failed.\n");
        return 1;
    }
    R.m = n;
    R.n = n;

    Tmp.m = n;
    Tmp.n = n;
    if (!(Tmp.mat = (DoubleComplex *)
                     malloc (n * n * sizeof (DoubleComplex)))) {
        printf ("brdErrChk: Malloc for Tmp.mat failed.\n");
        return 1;
    }
    /* R = Q' * A * (P')'
     * Q, A and P are n-by-n matrics, so call bmmult () not mmult ()
     */
    bmmult ("+", &Q, "H", &A, "N", NULL, &Tmp);
    bmmult ("+", &Tmp, "N", &P, "H", NULL, &R);

    mat = R.mat;
    for (j = 0; j < n; j++) {
        for (i = 0; i < n; i++) {
            if (i == j) {
                /* diagonal elements of B.      */
                mat[j * n + i].r -= d.vec[i];
            }
            else if (i == j - 1) {
                /* superdiagonal elements of B. */
                mat[j * n + i].r -= b.vec[i];
            }
        }
    }

    err = matnorm (&R) / (double)(n * n);
    printf ("Error in factorization: %1.3e\n", err);

    free (R.mat);
    free (Tmp.mat);
    return 0;
}