/****************************************************************************
 *
 * BlkLanApp.c
 *
 * Abstract: 
 *   This is the sample application of the block Lanczos tridiagonalization 
 *   of complex symmetric matrices. User can refer to this program to 
 *   invoke the routine to perform the tridiagonalization of a complex 
 *   symmetric matrix.
 *
 * Student:  Guohong Liu
 *
 * ID:       0385117
 *
 ****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <time.h>
#include "zcstrd.h"
#include "bidiagonal.h"
#include "blklanapp.h"

#ifdef MATLAB_COMP
doublecomplex *RandMat = NULL; 
doublecomplex *randVec = NULL; 
#endif

/***************************************************************************
 *
 * Abstract: 
 *   Main program.
 *   (1) Compiled by MATLAB_COMP, then command line to execute:
 *          BlkLanApp datafile. 
 *       where datafile is file contained data edited from MATLAB.
 *   (2) Compiled by INFO, or no compile preprocessor, then command line:
 *          BlkLanApp number1 number2
 *       where number1 is the order of matrix A, number2 is block size.
 *
 * Input   
 *   argc   Number of arguments.
 *   argv   Command line arguments.
 *
 * Output  
 *   None
 *
 * Return
 *   0      Normal exit
 *   1      Excetion was thrown
 *
 ***************************************************************************/

int main (int argc, char *argv[])
{
    /* Data needed by block Lanczos tridiagonalization routine. */
    int n;
    int bs;
    doublecomplex *A = NULL;
    doublecomplex *S = NULL;
    doublecomplex *Q = NULL;
    doublecomplex *P = NULL;
    doublecomplex *a = NULL;
    doublereal *b = NULL;
    doublereal *d = NULL;
    int info;
    char separator[20] = "===================\0";

    clock_t start, finish;
    float blkLanTime, brdTime;
#ifdef MATLAB_COMP
    DataFile dataFile;

    if (argc != 2) {
        printf ("Sample usage: BlkLanApp data.txt\n");
        printf ("data.txt: file which contains matrix and vector.\n");
        return 1;
    }
#else
    if (argc != 3) {
        printf ("Sample usage: BlkLanApp 256 4\n");
        printf ("256: size of matrix. 4: block size.\n");
        return 1;
    }
    if (atoi(argv[1]) % atoi(argv[2]) != 0) {
        printf ("Size of matrix should be dividable by block size.\n");
        return 1;
    }
#endif

    /* Allocate memories for block Lanczos tridiagonalization routine.  */
    if (blkLanInit (argv, &n, &bs, &A, &S, &Q, &P, &a, &b)) {
        return 1;
    }
    printf ("\n%s Block Lanczos Algorithm %s\n", separator, separator);
    start = clock();

    /* Perform block Lanczos routine.   */
    zcstrd_ (&n, A, &bs, S, a, b, Q, P, &info);
    if (info != 0) {
       return 1;
    }    
    
    finish = clock();
    blkLanTime = (double)(finish - start) / CLOCKS_PER_SEC;
    printf ("\nBlock Lanczos tridiagonalization: %5.2f seconds\n", 
            blkLanTime); 
    /* Check errors.                    */
    triErrChk (n, A, a, b, Q, P);
    freeBlkLan (S, Q, P, a, b);

    printf ("\n%s LAPACK Bidiagonalization %s\n", separator, separator);
    /* Allocate memories for bidiagonalization routine. */
    if (brdInit (n, A, &Q, &P, &d, &b)) {
        return 1;
    }
    start = clock();
    
    /* Perform bidiagonalization.       */
    brd (n, A, d, b, Q, P);
    
    finish = clock();
    brdTime = (double)(finish - start) / CLOCKS_PER_SEC;
    printf ("\nBidiagonalization: %5.2f seconds\n", brdTime);
    if (max(blkLanTime, brdTime) > 1.0) {
        printf (
        "\nRun time comparison: Block Lanczos/Bidiagonalization = %5.1f%%\n", 
        (blkLanTime / max(blkLanTime, brdTime)) * 100);
    }
    /* Check errors.                    */
    brdErrChk (n, A, d, b, Q, P);
    freeBrd (A, Q, P, d, b);

    return 0;
}

/***************************************************************************
 *
 * Abstract: 
 *   Allocate memories for block Lanczos algorithm routine.
 *
 * Input
 *   argv       Parameters of command line
 *
 * Output
 *   order      Order of matrix A.
 *   blkSize    Block size.
 *   A          Pointer to the memory of matrix A.
 *   S          Pointer to the starting matrix S.
 *   Q          Pointer to the memory of matrix Q.
 *   P          Pointer to the memory of matrix P.
 *   a          Pointer to the memory of vector a.
 *   b          Pointer to the memory of vector b.
 *
 * Return
 *   0      Successful exit
 *   1      Memory allocation exception occurs
 *
 ***************************************************************************/

int blkLanInit (char *argv[], int *order, int *blkSize, 
                doublecomplex **A, doublecomplex **S, 
                doublecomplex **Q, doublecomplex **P, 
                doublecomplex **a, doublereal **b)
{ 
#ifdef MATLAB_COMP
    DataFile dataFile;
#endif
    int n, bs;
 
#ifdef MATLAB_COMP
    dataFile = readDataFile (argv[1]);
    n = dataFile.n;
    bs = dataFile.bs;
    *order = n;
    *blkSize = bs;
    *A = dataFile.A;
    *S = dataFile.S;
    RandMat = dataFile.RandMat;
    randVec = dataFile.randVec;
#else
    n = atoi (argv[1]);
    bs = atoi (argv[2]);
    *order = n;
    *blkSize = bs;
    if (!(*A = (doublecomplex *)
                malloc (n * n * sizeof (doublecomplex)))){
        printf ("blkLanInit: Malloc for *A failed.\n");
        return 1;
    }
    /* Before generating random number, it is necessory to initialize
     * the random seed first.
     */
    initRandSeed ();
    randComplexSymMat (n, *A);

    if (!(*S = (doublecomplex *)
                malloc (n * bs * sizeof (doublecomplex)))){
        blkLanInitExcept (*A, *S, *Q, *P, *a, *b);
        printf ("blkLanInit: Malloc for *S failed.\n");
        return 1;
    }
    randOrthMat (n, bs, *S);
#endif

    if (!(*Q = (doublecomplex *)
                malloc (n * n * sizeof (doublecomplex)))) {
        blkLanInitExcept (*A, *S, *Q, *P, *a, *b);
        printf ("blkLanInit: Malloc for *Q failed.\n");
        return 1;
    }
    if (!(*P = (doublecomplex *)
                malloc (n * n * sizeof (doublecomplex)))) {
        blkLanInitExcept (*A, *S, *Q, *P, *a, *b);
        printf ("blkLanInit: Malloc for *P failed.\n");
        return 1;
    }

    if (!(*a = (doublecomplex *)malloc (n * sizeof (doublecomplex)))) {
        blkLanInitExcept (*A, *S, *Q, *P, *a, *b);
        printf ("blkLanInit: Malloc for *a failed.\n");
        return 1;
    }
    if (!(*b = (doublereal *)malloc ((n - 1) * sizeof (doublereal)))) {
        blkLanInitExcept (*A, *S, *Q, *P, *a, *b);
        printf ("blkLanInit: Malloc for *b failed.\n");
        return 1;
    }
    return 0;
}

/***************************************************************************
 *
 * Abstract
 *   Initialize for the bidiagonalization of a complex symmetric matrix.  
 *
 * Input
 *   n      Order of matrix A.
 *
 * Output
 *   A      Pointer to the memory of matrix A.
 *   Q      Pointer to the memory of matrix Q.
 *   P      Pointer to the memory of matrix P.
 *   d      Pointer to the memory of vector d.
 *   b      Pointer to the memory of vector b.
 *
 * Return
 *   0      Successful exit
 *   1      Memory allocation exception occurs
 *
 ***************************************************************************/

int brdInit (int n, doublecomplex *A, doublecomplex **Q, doublecomplex **P, 
             doublereal **d, doublereal **b)
{
    /* Allocate memory for the unitary matrices Q, and P.    */
    if (!(*Q = (doublecomplex *)malloc (n * n * sizeof (doublecomplex)))) {
        printf ("brdInit: Malloc for *Q failed.\n");
        return 1;
    }
    if (!(*P = (doublecomplex *)malloc (n * n * sizeof (doublecomplex)))) {
        brdInitExcept (A, *Q, *P, *d, *b);
        printf ("brdInit: Malloc for *P failed.\n");
        return 1;
    }

    /* Allocate memory for diagonal elements.               */
    if (!(*d = (doublereal *)malloc (n * sizeof (doublereal)))) {
        brdInitExcept (A, *Q, *P, *d, *b);
        printf ("brdInit: Malloc for *d failed.\n");
        return 1;
    }

    /* Allocate memory for off-diagonal elements.           */
    if (!(*b = (doublereal *)malloc ((n - 1) * sizeof (doublereal)))) {
        brdInitExcept (A, *Q, *P, *d, *b);
        printf ("brdInit: Malloc for *b failed.\n");
        return 1;
    }

    return 0;
}

/***************************************************************************
 *
 * Abstract: 
 *   Handle exceptions thrown in the intialization of block Lanczos
 *   tridiagonalization. Release memories one by one.
 *
 * Input   
 *   A          Pointer to the memory of matrix A.
 *   S          Pointer to the starting matrix S.
 *   Q          Pointer to the memory of matrix Q.
 *   P          Pointer to the memory of matrix P.
 *   a          Pointer to the memory of vector a.
 *   b          Pointer to the memory of vector b.
 *
 * Output
 *   None
 *
 ***************************************************************************/

void blkLanInitExcept (doublecomplex *A, doublecomplex *S, 
                       doublecomplex *Q, doublecomplex *P, 
                       doublecomplex *a, doublereal *b)
{
    if (A != NULL) {
        free (A);
        A = NULL;
    }
    if (S != NULL) {
        free (S);
        S = NULL;
    }
    if (Q != NULL) {
        free (Q);
        Q = NULL;
    }
    if (P != NULL) {
        free (P);
        P = NULL;
    }
    if (a != NULL) {
        free (a);
        a = NULL;
    }
    if (b != NULL) {
        free (b);
        b = NULL;
    }
}

/***************************************************************************
 *
 * Abstract: 
 *   Handle exceptions thrown in the intialization of bidiagonalization.
 *   Release memories one by one.
 *
 * Input   
 *   A          Pointer to the memory of matrix A.
 *   S          Pointer to the starting matrix S.
 *   Q          Pointer to the memory of matrix Q.
 *   P          Pointer to the memory of matrix P.
 *   a          Pointer to the memory of vector a.
 *   b          Pointer to the memory of vector b.
 *
 * Output
 *   None
 *
 ***************************************************************************/

void brdInitExcept (doublecomplex *A, doublecomplex *Q, doublecomplex *P, 
                    doublereal *d, doublereal *b)
{
    if (A != NULL) {
        free (A);
        A = NULL;
    }
    if (Q != NULL) {
        free (Q);
        Q = NULL;
    }
    if (P != NULL) {
        free (P);
        P = NULL;
    }
    if (d != NULL) {
        free (d);
        d = NULL;
    }
    if (b != NULL) {
        free (b);
        b = NULL;
    }
}

/***************************************************************************
 *
 * Abstract
 *   Free memories allocated for block Lanczos tridiagonalization.
 *
 * Input   
 *   S          Pointer to the starting matrix S.
 *   Q          Pointer to the memory of matrix Q.
 *   P          Pointer to the memory of matrix P.
 *   a          Pointer to the memory of vector a.
 *   b          Pointer to the memory of vector b.
 *
 * Output
 *   None
 *
 ***************************************************************************/

void freeBlkLan (doublecomplex *S, doublecomplex *Q, doublecomplex *P, 
                 doublecomplex *a, doublereal *b)
{
    /* Free memories except matrix A because it is used later 
     * for bidiagonalization comparison.
     */
    free (S);
    S = NULL;
    free (Q);
    Q = NULL;
    free (P);
    P = NULL;
    free (a);  
    a = NULL;
    free (b);
    b = NULL;
#ifdef MATLAB_COMP
    free (RandMat);
    free (randVec);
#endif
}

/***************************************************************************
 *
 * Abstract
 *   Free memories allocated for bidiagonalization.
 *
 * Input   
 *   A          Pointer to the memory of matrix A.
 *   Q          Pointer to the memory of matrix Q.
 *   P          Pointer to the memory of matrix P.
 *   a          Pointer to the memory of vector a.
 *   b          Pointer to the memory of vector b.
 *
 * Output
 *   None
 *
 ***************************************************************************/

void freeBrd (doublecomplex *A, doublecomplex *Q, doublecomplex *P, 
              doublereal *d, doublereal *b)
{
    /* Release all memories.  */
    free (A);
    A = NULL;
    free (Q);
    Q = NULL;
    free (P);
    P = NULL;
    free (d);
    d = NULL;
    free (b);
    b = NULL;
}