/****************************************************************************
 *
 * LapackWrap.c
 *
 * Abstract 
 *   Implement routines of LAPACK and BLAS wrapper.
 *
 * Student:  Guohong Liu
 *
 * ID:       0385117
 *
 ****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <time.h>

#include "f2c.h"
#include "blaswrap.h"
#include "fblaswr.h"
#include "computing.h"

/***************************************************************************
 *
 * Global variablies
 *
 *   For the sake of performance, specify work buffers which frequently
 *   involve in computing:
 *   *mat     Work buffer for a matrix. Max size = n * bs 
 *   *tau     Work buffer for QR fact.  Max size = n 
 *   *work    Work buffer for QR fact.  Max size = n
 * 
 ***************************************************************************/

DoubleComplex *mat = NULL;     
DoubleComplex *tau = NULL;    
DoubleComplex *work = NULL;   

integer iseed[4];
/*  Flag to specify that random number is uniform (-1,1)    */
int idist = 2;

/***************************************************************************
 *
 * Abstract 
 *   Allocate work memory for routines of Lapack wrapper. It is necessory
 *   to intialize the Lapack wrapper before calling routines of this file,
 *   and free the wrapper when computation is done.
 *   
 * Input    
 *   n      The order of matrix involved in the compuation of wrapper.
 *   bs     The block size.
 *
 * Output     
 *   None
 * 
 * Return
 *   0      Successful exit
 *   1      Memory allocation exception occurs
 *
 ***************************************************************************/

int wrapperInit (int n, int bs)
{
    if (!(mat = (DoubleComplex *)
                malloc (n * bs * sizeof (DoubleComplex)))) {
        printf ("wrapperInit: Malloc for *mat failed.\n");
        return 1;
    }
    if (!(tau = (DoubleComplex *)malloc (n * sizeof (DoubleComplex)))) {
        printf ("wrapperInit: Malloc for *tau failed.\n");
        free (mat);
        mat = NULL;
        return 1;
    }
    if (!(work = (DoubleComplex *)malloc (n * sizeof (DoubleComplex)))) {
        free (mat);
        mat = NULL;
        free (tau);
        tau = NULL;
        printf ("wrapperInit: Malloc for *work failed.\n");
        return 1;
    }
    
    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Release work memories occupied by LAPACK and BLAS wrapper.
 *
 * Input    
 *   None
 *
 * Output
 *   None
 *
 ***************************************************************************/

void freeWrapper ()
{
     free (mat);
     mat = NULL;
     free (tau);
     tau = NULL;
     free (work);
     work = NULL;
}

/***************************************************************************
 *
 * Abstract 
 *   QR factorization. A = Q * R.
 *  
 * Input   
 *   A      Matrix. Size = m * n.    
 *          
 * Output  
 *   Q      Matrix. Size = m * n.
 *   R      Matrix. Size = n * n.
 *
 * Return  
 *   0      Success exit
 *   1      Memory allocation exception occurs.
 *
 * NOTE
 *   For the sake of performance, this routine does the QR 
 *   factorization for only n-by-bs matrices, where n is the order
 *   of matrix A which is given the block Lanczos routine to reduce,
 *   and bs is the block size. This routine is not suitable to 
 *   do QR factorization for n-by-n square matrices.
 *
 *
 ***************************************************************************/

int qr (DoubleComplexMat *A, DoubleComplexMat *Q, DoubleComplexMat *R)
{
    int m = A->m;
    int n = A->n;
    int dim = n;
    int lwork = n;  
    int info;
    DoubleComplex alpha, beta;

    /* Matrix A will be overwriten in QR, so make a copy of it.   */
    zlacpy_ ("", &m, &n, A->mat, &m, mat, &m);

    /* R is lower triangular. Initialize R so that its elements in 
     * lower triangular are all 0.  
     */
    alpha.r = 0;
    alpha.i = 0;
    beta.r = 0;
    beta.i = 0;
    zlaset_ ("L", &n, &n, &alpha, &beta, R->mat, &n);
    
    zgeqrf_ (&m, &n, mat, &m, tau, work, &lwork, &info);
    zlacpy_ ("U", &n, &n, mat, &m, R->mat, &n);

    zungqr_ (&m, &n, &dim, mat, &m, tau, work, &lwork, &info);
    zlacpy_ ("", &m, &n, mat, &m, Q->mat, &m);

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Apply Gram-Schmidt method to orthogonalizing A against B.
 *
 * Input    
 *   A      Matrix A
 *   B      Matrix B
 *
 * Output     
 *   A      The updated matrix
 *
 ***************************************************************************/

void orthogonalize (DoubleComplexMat *A, DoubleComplexMat *B)
{
    DoubleComplex *v_A, *v_B, ret;
    int i, j, cola, colb, n, inc = 1;
    
    n = A->m;
    /* Set v_A as the first column (vector) of the matrix A.        */
    v_A = A->mat;
    for (i = 1; i <= A->n; i++) {
        /* Set v_B as the first column (vector) of the matrix B.    */
        v_B = B->mat;
        
        for (j = 1; j <= B->n; j++) {
            /* Compute A(:,i) = A(:,i) - (B(:,j)'*A(:,i)) * B(:,j)  
             * First compute ret = B(:,j)' * A(:,i) 
             */
            zdotc_ (&ret, &n, v_B, &inc, v_A, &inc);
            ret.r = -ret.r;
            ret.i = -ret.i;
            zaxpy_ (&n, &ret, v_B, &inc, v_A, &inc);

            /* Set v_B as the next column of matrix B.  */
            v_B += B->m;
        }

        /* Set v_A as the next column of matrix A.      */
        v_A += A->m;
    }
}

/***************************************************************************
 *
 * Abstract
 *   Compute the multiplication and addition of matrices, 
 *      R =  op(A) * op(B) + C,  
 *      R = -op(A) * op(B) + C, 
 *   BLAS doesn't support conjg(B), we must implement it.
 *
 * Input   
 *   signa   "+", op(A) * op(B) + C
 *           "-", -op(A) * op(B) + C
 *   A       Matrix. Size of op(A) = m * k          
 *   transa  "N", op(A) = A 
 *           "T", op(A) = A.'          nonconjugate transpose.  
 *           "H", op(A) = A'           conjugate transpose  
 *   B       Matrix. Size of op(B) = k * n
 *   transb  "N", op(B) = B 
 *           "T", op(B) = B.'          nonconjugate transpose.  
 *           "H", op(B) = B'           conjugate transpose  
 *           "C", op(B) = conjg (B)    just conjugate. 
 *   C       Matrix. Size = m * n
 *           if C = NULL, then R = (+, -) op(A) * op(B)
 *
 * Output  
 *   R       Matrix. Size = m * n.    
 *
 * Return
 *   0       Successful exit
 *   1       Exception occurs
 * 
 * NOTE
 *   (1) The function assumes matrix R is not matrix A or B to 
 *       achieve the best performance. That means the result of 
 *       computing should be put in C or another matrix. 
 *   (2) The function doesn't deal with the case that op(A) = "C" 
 *       because block Lanczos method doesn't has this operation.
 *   (3) Only one of the matrices, A, B or C, can be n-by-n size 
 *       if one of them involves "T", "H" or "C" operation. 
 *       Otherwise call routine bmmult (), where 'b' means big 
 *       matrix-matrix operation. If there is no such operations 
 *       at all, matrix A and B  can be n-by-n size.
 * 
 ***************************************************************************/

int mmult (char *signa, DoubleComplexMat *A, char *transa, 
           DoubleComplexMat *B, char *transb, 
           DoubleComplexMat *C, DoubleComplexMat *R)
{
    DoubleComplex alpha, beta;
    DoubleComplex *B_mat;
    char opa, opb; 
    int m, n, k;
    int incx = 1, len;

    if (signa[0] == '+') {
        alpha.r = 1;
        alpha.i = 0;
    }
    else if (signa[0] == '-') {
        alpha.r = -1;
        alpha.i = 0;
    }
    else {
        printf ("mmult: Wrong sign before matrix A.\n");
        return 1;
    }

    if (C == NULL) {
        beta.r = 0;
        beta.i = 0;
    }
    else {
        beta.r = 1;
        beta.i = 0;
    }
    
    if (transb[0] == 'C') {
        B_mat = mat;
        /* Copy matrix B to memory B_mat because zlacgv_ () puts 
         * the solution in B_mat 
         */
        zlacpy_ ("", &(B->m), &(B->n), B->mat, &(B->m), B_mat, &(B->m));
        len = B->m * B->n;
        /* Compute conjugate of B   */
        zlacgv_ (&len, B_mat, &incx);
    }
    else {
        B_mat = B->mat;
    }

    if ((C != NULL) && (R->mat != C->mat)) {
        /* Caller specifies that the result will not be stored in
         * C, but in R. Copy matrix C to R.
         */
        zlacpy_ ("", &(C->m), &(C->n), C->mat, &(C->m), R->mat, &(R->m));
    }

    /* Set m, n, k, opa, opb        */
    if ((transa[0] == 'N') || (transa[0] == 'C')) {
        opa = 'N';
        m = A->m;
        k = A->n;
    }
    else if (transa[0] == 'T') {
        opa = 'T';
        m = A->n;
        k = A->m;
    }
    else if (transa[0] == 'H') {
        opa = 'C';
        m = A->n;
        k = A->m;
    }

    if ((transb[0] == 'N') || (transb[0] == 'C')) {
        opb = 'N';
        n = B->n;
    }
    else if (transb[0] == 'T') {
        opb = 'T';
        n = B->m;
    }
    else if (transb[0] == 'H') {
        opb = 'C';
        n = B->m;
    }
   
    zgemm_ (&opa, &opb, &m, &n, &k, &alpha, A->mat, &(A->m),
            B_mat, &(B->m), &beta, R->mat, &(R->m));

    return 0;
}

/***************************************************************************
 *
 * Abstract
 *   Compute multiplication and addition of matrices, 
 *      R =  op(A) * op(B) + C,  
 *      R = -op(A) * op(B) + C, 
 *   BLAS doesn't support op(B) = conjg (B).
 * 
 * Input   
 *   signa   "+", op(A) * op(B) + C
 *           "-", -op(A) * op(B) + C
 *   A       Matrix. Size of op(A) = m * k          
 *   transa  "N", op(A) = A 
 *           "T", op(A) = A.'          nonconjugate transpose.  
 *           "H", op(A) = A'           conjugate transpose  
 *   B       Matrix. Size of op(B) = k * n
 *   transb  "N", op(B) = B 
 *           "T", op(B) = B.'          nonconjugate transpose.  
 *           "H", op(B) = B'           conjugate transpose  
 *           "C", op(B) = conjg (B)    just conjugate.
 *   C       Matrix. Size = m * n
 *           if C = NULL, then R = (+, -) op(A) * op(B)
 *
 * Output  
 *   R       Matrix. Size = m * n.   
 *
 * Return
 *   0       Successful exit
 *   1       Exception occurs
 *
 * NOTE
 *   This function is the same as mmult () except matrices A and B are two 
 *   bigger matrices, n-by-n for example. mmult () has better performance
 *   than this one.
 *          
 ***************************************************************************/

int bmmult (char *signa, DoubleComplexMat *A, char *transa, 
            DoubleComplexMat *B, char *transb, 
            DoubleComplexMat *C, DoubleComplexMat *R)
{
    DoubleComplex alpha, beta;
    DoubleComplexMat B1;
    char opa, opb; 
    int m, n, k;

    if (signa[0] == '+') {
        alpha.r = 1;
        alpha.i = 0;
    }
    else if (signa[0] == '-') {
        alpha.r = -1;
        alpha.i = 0;
    }
    else {
        printf ("mmult: Wrong sign before matrix A.\n");
        return 1;
    }

    if (C == NULL) {
        beta.r = 0;
        beta.i = 0;
    }
    else {
        beta.r = 1;
        beta.i = 0;
    }
    
    if (transb[0] == 'C') {
        /* Compute conjugate of B   */
        if (!(B1.mat = (DoubleComplex *)
                        malloc (B->m * B->n * sizeof (DoubleComplex)))) {
            printf ("bmmult: Malloc for B1.mat failed.\n");
            return 1;
        }
        B1.m = B->m;
        B1.n = B->n;
        conjmat (B, &B1);
    }
    else {
        B1.m = B->m;
        B1.n = B->n;
        B1.mat = B->mat;
    }

    if ((C != NULL) && (R->mat != C->mat)) {
        /* Caller specifies that the result will not be stored in
         * C, but in R. Copy matrix C to R.
         */
        zlacpy_ ("", &(C->m), &(C->n), C->mat, &(C->m), R->mat, &(R->m));
    }

    /* Set m, n, k, opa, opb        */
    if ((transa[0] == 'N') || (transa[0] == 'C')) {
        opa = 'N';
        m = A->m;
        k = A->n;
    }
    else if (transa[0] == 'T') {
        opa = 'T';
        m = A->n;
        k = A->m;
    }
    else if (transa[0] == 'H') {
        opa = 'C';
        m = A->n;
        k = A->m;
    }

    if ((transb[0] == 'N') || (transb[0] == 'C')) {
        opb = 'N';
        n = B1.n;
    }
    else if (transb[0] == 'T') {
        opb = 'T';
        n = B1.m;
    }
    else if (transb[0] == 'H') {
        opb = 'C';
        n = B1.m;
    }
   
    zgemm_ (&opa, &opb, &m, &n, &k, &alpha, A->mat, &(A->m),
             B1.mat, &(B1.m), &beta, R->mat, &(R->m));

    if (transb[0] == 'C') {
        free (B1.mat);
    }

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Square matrix-matrix left division. R = A / B, where B is 
 *   upper triangle. 
 * 
 * Input    
 *   A      Matrix. Size = n * n  
 *   B      Upper triangle matrix. Size = n * n
 *
 * Output   
 *   R      Matrix. Size = n * n.  
 *
 * Return  
 *   0      Successful exit
 *   1      Exception occurs
 * 
 * NOTE
 *   The function doesn't deal with the case that B = A / B, ie. the 
 *   result is stored in matrix B.
 *
 ***************************************************************************/

int mdiv (DoubleComplexMat *A, DoubleComplexMat *B, DoubleComplexMat *R)
{
    int n = A->m;
    DoubleComplex alpha;
    int info;

    if (R->mat != A->mat) {
        /* Copy matrix A to matrix Result because ztrsm_ () 
         * puts the solution in A 
         */
        zlacpy_ ("", &n, &n, A->mat, &n, R->mat, &n);
    }

    alpha.r = 1;
    alpha.i = 0;
    ztrsm_ ("R", "U", "N", "N", &n, &n, &alpha, B->mat, &n, R->mat, &n);

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute addition of two matrices,
 *      R = A + B       
 *   
 * Input    
 *   A      Matrix. Size = m * n.
 *   B      Matrix. Size = m * n.
 *
 * Output   
 *   R      Matrix. Size = m * n.   
 *
 * Return  
 *   0      Successful exit
 *   1      Exception occurs
 * 
 * NOTE
 *   The function is OK for the case B = A + B, but it doesn't 
 *   deal with the case that A = A + B, ie. the result is stored in A. 
 *
 ***************************************************************************/

int mplus (DoubleComplexMat *A, DoubleComplexMat *B, DoubleComplexMat *R)
{
    DoubleComplex alpha;
    integer n = (A->m) * (A->n);
    integer inca = 1, incx = 1;

    if (R->mat != B->mat) {
        /* Copy matrix B to matrix Result because zaxpy_ () 
         * puts the solution in B 
         */
        zlacpy_ ("", &(B->m), &(B->n), B->mat, &(B->m), R->mat, &(R->m));
    }

    alpha.r = 1;
    alpha.i = 0;
    /* Because matrix is stored as one dimension array. Treat the 
     * matrix as vector for matrix-matrix addition.     
     */
    zaxpy_ (&n, &alpha, A->mat, &inca, R->mat, &incx);

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute the conjugate of a matrix, R = A.'
 *
 * Input    
 *   A      Matrix A
 *
 * Output     
 *   R      Matrix R
 *
 * NOTE
 *   This function assumes the result matrix R is different from
 *   matrix A for the sake of performance. The result is not stored 
 *   in previous matrix. 
 *
 ***************************************************************************/

void conjmat (DoubleComplexMat *A, DoubleComplexMat *R)
{

    int n = (A->m) * (A->n);
    int incx = 1;

    /* Copy matrix A to R because zlacgv_ () puts the solution in A */
    zlacpy_ ("", &(A->m), &(A->n), A->mat, &(A->m), R->mat, &(R->m));
    zlacgv_ (&n, R->mat, &incx);
}

/***************************************************************************
 *
 * Abstract 
 *   Scale a matrix. It computes
 *      R = alpha * A,  where alpha is a real constant. 
 *   
 * Input    
 *   alpha      double number alpha
 *   A          Matrix A
 *
 * Output     
 *   R          Matrix R
 *
 * NOTE
 *   The function assumes the result is not stored in matrix A,
 *   but in a second matrix, R. 
 *
 ***************************************************************************/

void mscal (double alpha, DoubleComplexMat *A, DoubleComplexMat *R)
{
    integer n = (A->m) * (A->n);
    integer inc = 1;

    /* Copy matrix A to R, use R in BLAS function zaxpy_ (). */
    zlacpy_ ("", &(A->m), &(A->n), A->mat, &(A->m), R->mat, &(R->m));

    /* Because matrix is stored as one dimension array. Treat the 
     * matrix as vector for matrix-matrix addition.     
     */
    zdscal_ (&n, &alpha, R->mat, &inc);
}

/***************************************************************************
 *
 * Abstract 
 *   Compute the norm of a matrix A
 *
 * Input    
 *   A      Matrix A
 *      
 * Output   
 *   None
 * 
 * Return
 *   norm of matrix A
 *
 ***************************************************************************/

double matnorm (DoubleComplexMat *A)
{
    double err = zlange_ ("F", &(A->m), &(A->n), A->mat, &(A->m), NULL);
    
    return err;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute multiplication of a matrix and a vector,
 *      r = op(A) * x + y
 *   where A is not a complex symmetric matrix.
 * 
 * Input    
 *   A       Complex matrix. Size = m * n
 *   transa  "N", op(A) = A 
 *           "T", op(A) = A.'  nonconjugate transpose.  
 *   x       Vector. Size = n
 *   y       Vector. Size = n
 *      
 * Output     
 *   r       Vector. Size = n
 *
 * NOTE
 *   The function assumes the result is stored in vector r which
 *   is different from vector x. That is OK to store the result in 
 *   vector y.
 *  
 ***************************************************************************/

int mvmult (DoubleComplexMat *A, char *transa, DoubleComplexVec *x, 
            DoubleComplexVec *y, DoubleComplexVec *r)
{
    DoubleComplex alpha, beta;
    int inc = 1, columns = 1;

    if (y != NULL) {
        if (r->vec != y->vec) {
            /* Copy vector y to vector result.  */
            zlacpy_ ("", &(y->n), &columns, y->vec, &(y->n), 
                     r->vec, &(r->n));
        }
        beta.r = 1;
        beta.i = 0;
    }
    else {
        beta.r = 0;
        beta.i = 0;
    }

    alpha.r = 1;
    alpha.i = 0;
    
    zgemv_ (transa, &(A->m), &(A->n), &alpha, A->mat, &(A->m), 
            x->vec, &inc, &beta, r->vec, &inc);  

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute multiplication of a matrix and a vector,
 *      r = A * x + y
 *   where A is a complex symmetric matrix.
 *
 * Input    
 *   A       Complex symmetric matrix. Size = m * n
 *   x       Vector. Size = n
 *   y       Vector. Size = n
 *      
 * Output     
 *   r       Vector. Size = n
 *
 * NOTE
 *   The function assumes the result is stored in vector r which
 *   is different from vector x. That is OK to store the result in 
 *   vector y.
 *
 ***************************************************************************/

int symvmult (DoubleComplexMat *A, DoubleComplexVec *x, 
              DoubleComplexVec *y, DoubleComplexVec *r)
{
    DoubleComplex alpha, beta;
    int inc = 1, columns = 1;

    if (y != NULL) {
        if (r->vec != y->vec) {
            /* Copy vector y to vector result.  */
            zlacpy_ ("", &(y->n), &columns, y->vec, &(y->n), 
                        r->vec, &(r->n));
        }
        beta.r = 1;
        beta.i = 0;
    }
    else {
        beta.r = 0;
        beta.i = 0;
    }

    alpha.r = 1;
    alpha.i = 0;
    
    zsymv_ ("U", &(A->m), &alpha, A->mat, &(A->m), 
            x->vec, &inc, &beta, r->vec, &inc);  

    return 0;
}

/***************************************************************************
 *
 * Abstract    
 *   Compute the multiplication of two vectors,
 *      r = x.' * y
 *      r = conjg(x.') * y
 *
 * Input   
 *   x       Vector.
 *   transx  "T", op(x) = x.'          nonconjugate transpose.  
 *           "H", op(a) = conjg (x.'). conjugate transpose  
 *   y       Vector.
 *
 * Output
 *   r       Double complex number.
 *
 ***************************************************************************/

int vmult (DoubleComplexVec *x, char *transx, DoubleComplexVec *y, 
           DoubleComplex *r)
{

    int inc = 1;

    if (transx[0] == 'T') {
        /* Compute r = x.' * y  */
        zdotu_ (r, &(x->n), x->vec, &inc, y->vec, &inc);
    }
    else if (transx[0] == 'H') {
        /* Compute r = x' * y   */
        zdotc_ (r, &(x->n), x->vec, &inc, y->vec, &inc);
    }
    else {
        printf ("vmult: Wrong parameter for transx.\n");
        return 1;
    }

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Constant times a vector plus a vector, 
 *      y = alpha * x + y, 
 *      y = -alpha * x + y
 *   where alpha is a complex scalor.
 *
 * Input    
 *   alpha      Complex scalor
 *   x          Complex vector
 *   y          Complex vector
 *
 * Output     
 *
 * NOTE
 *   The function assumes the result is stored in vector y.
 *
 ***************************************************************************/

int vztplus (char *signx, DoubleComplex *alpha, DoubleComplexVec *x, 
             DoubleComplexVec *y, DoubleComplexVec *r)
{
    int columns = 1, inc = 1;
    DoubleComplex alpha1;

    if (signx[0] == '+') {
        alpha1.r = alpha->r;        
        alpha1.i = alpha->i;
    }
    else if (signx[0] == '-') {
        alpha1.r = -alpha->r;
        alpha1.i = -alpha->i;
    }
    else {
        printf ("vplus: Wrong sign before vector x.\n");
        return 1;
    }

    zaxpy_ (&(x->n), &alpha1, x->vec, &inc, r->vec, &inc);

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Constant times a vector plus a vector, 
 *      y = alpha * x + y, 
 *      y = -alpha * x + y
 *   where alpha is a real scalor.
 *
 * Input    
 *   alpha      Real scalor
 *   x          Complex vector
 *   y          Complex vector
 *
 * Output     
 *   r          Complex vector
 *
 * NOTE
 *   The function assumes the result is stored in vector y.
 *
 ***************************************************************************/

int vdtplus (char *signx, DoubleReal alpha, DoubleComplexVec *x, 
             DoubleComplexVec *y, DoubleComplexVec *r)
{
    DoubleComplex alpha_c;
    int columns = 1, inc = 1;

    if (signx[0] == '+') {
        alpha_c.r = alpha;
        alpha_c.i = 0;
    }
    else if (signx[0] == '-') {
        alpha_c.r = -alpha;
        alpha_c.i = 0;
    }
    else if (signx[0] != '+') {
        printf ("vplus: Wrong sign before vector x.\n");
        return 1;
    }

    zaxpy_ (&(x->n), &alpha_c, x->vec, &inc, r->vec, &inc);

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   A real vector plus a real vector, y = x + y, 
 *
 * Input    
 *   x          Real vector
 *   y          Real vector
 *
 * Output   
 *   r          Real vector
 *
 * NOTE
 *   The function assumes the result is stored in vector y.
 *   
 ***************************************************************************/

int vplus (DoubleRealVec *x, DoubleRealVec *y, DoubleRealVec *r)
{
    DoubleReal alpha = 1.0;
    int columns = 1, inc = 1;

    daxpy_ (&(x->n), &alpha, x->vec, &inc, r->vec, &inc);

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Scale the vector. It computes,
 *      r = alpha * x,   
 *   where alpha is a real scalor.
 *
 * Input    
 *   alpha      double number alpha
 *   x          vector x.
 *
 * Output     
 *   r          vector r
 *
 * NOTE
 *   The function assumes the result is stored in vector r which
 *   is different from vector x. 
 *
 ***************************************************************************/

void vscal (double alpha, DoubleComplexVec *x, DoubleComplexVec *r)
{
    int inc = 1, columns = 1;

    /* Copy vector x to r, use r in BLAS function zaxpy_ (). */
    zlacpy_ ("", &(x->n), &columns, x->vec, &(x->n), r->vec, &(r->n));
    zdscal_ (&(x->n), &alpha, r->vec, &inc);
}

/***************************************************************************
 *
 * Abstract 
 *   Compute the norm of a vector.
 *
 * Input
 *   x      Vector x
 *
 * Output     
 *   None
 * 
 * Return
 *   norm of the vector x
 *
 ***************************************************************************/

double vecnorm (DoubleComplexVec *x)
{
    int columns = 1;
    double err = zlange_ ("F", &(x->n), &columns, x->vec, &(x->n), NULL);
    
    return err;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute the conjugate of one vector, r = x.'
 *
 * Input    
 *   x      Vector x
 *
 * Output     
 *   r      Vector r
 *
 * NOTE
 *   The function assumes the vector r is different from vector x. 
 * 
 ***************************************************************************/

void conjvec (DoubleComplexVec *x, DoubleComplexVec *r)
{
    int incx = 1, columns = 1;

    /* Copy vector x to r because zlacgv_ () puts the solution in x */
    zlacpy_ ("", &(x->n), &columns, x->vec, &(x->n), r->vec, &(r->n));
    zlacgv_ (&(r->n), r->vec, &incx);
}

/***************************************************************************
 *
 * Abstract 
 *   Get the epsilon of the machine.
 *
 * Input    
 *   None
 *
 * Output
 *   None
 *
 ***************************************************************************/

double getEps ()
{
    return dlamch_ ("Epsilon");
}

/***************************************************************************
 *
 * Abstract 
 *   Initialize the seed of the random number generator DLARNV.
 *   The seed array elements must be between 0 and 4095, otherwise 
 *   they will be reduced mod 4096, and iseed(4) must be odd. After 
 *   function zlarnv_ () returns a vector of n random complex numbers 
 *   from a uniform, the seed is updated. 
 *
 * Input    
 *   None
 *
 * Output   
 *   None        
 *
 ***************************************************************************/

void initRandSeed ()
{
    int i;
    
    for (i = 0; i < 3; i++) {
        /* Seed the random-number generator with current time so that
         * the numbers will be different every time we run.
         */
	    iseed[i] = (unsigned)time (NULL) % 4096;
    }
    iseed[3] = 1;
}

/***************************************************************************
 *
 * Abstract 
 *   Generate a random complex matrix. The real and imaginary part of 
 *   each element are in uniform (-1,1) distribution.  
 *
 * Input 
 *   A      Matrix A
 *
 * Output     
 *   A      Matrix A
 *
 ***************************************************************************/

void randComplexMat (DoubleComplexMat *A)
{
    int len;

    len = A->m * A->n;
    zlarnv_ (&idist, &iseed[0], &len, A->mat);
}
                  
/***************************************************************************
 *
 * Abstract 
 *   Generate a random complex symmetric matrix. The real and imaginary 
 *   part of each element are in uniform (-1,1) distribution.  
 * 
 * Input    
 *   n      the order of matrix A
 *   A_m    the memory to store matrix A
 *
 * Output     
 *   A_m    the random complex symmetric matrix.
 *
 ***************************************************************************/

void randComplexSymMat (int n, DoubleComplex *A_m)
{
    DoubleComplex *ptr;
    int len, i, j;

    /* Create the lower triangular part of the matrix.  */
    ptr = A_m;
    for (j = 0; j < n; j++) {
        len = n - j;
        zlarnv_ (&idist, &iseed[0] ,&len, ptr);
        ptr += len + j + 1;
    }

    /* Create the upper triangular part of matrix by 
     * copying the lower triangular part
     */
    ptr = A_m;
    for (j = 0; j < n; j++) {
        for (i = 0; i < j; i++) {
            ptr[j * n + i] = ptr[i * n + j];
        }
    }
}

/***************************************************************************
 *
 * Abstract 
 *   Generate starting matrix of orthonormal columns by using 
 *   QR factorizatoin.   
 *
 * Input
 *   m      the rows of matrix A
 *   n      the columns of matrix A
 *   A_m    the memory to store matrix A
 *
 * Output     
 *   A_m    the random complex symmetric matrix.
 *
 ***************************************************************************/

int randOrthMat (int m, int n, DoubleComplex *A_m)
{
    DoubleComplex *Rand;
    DoubleReal *RealMat;
    int len = m * n;
    int dim = min (m, n);
    DoubleComplex *tau, *work;
    int lwork = n;  
    int i, j, info;
    DoubleComplex alpha, beta;

    /* Generate a random real matrix. */
    if (!(RealMat = (DoubleReal *)malloc (m * n * sizeof (DoubleReal)))){
        printf ("randOrthComplexMat: Malloc for Rand.mat failed.\n");
        return 1;
    }
    dlarnv_ (&idist, &iseed[0] ,&len, RealMat);

    /* Copy the random real matrix to a complex matrix. */
    if (!(Rand = (DoubleComplex *)
                    malloc (m * n * sizeof (DoubleComplex)))){
        printf ("randOrthComplexMat: Malloc for Rand.mat failed.\n");
        return 1;
    }
    for (j = 0; j < n; j++) {
        for (i = 0; i < m; i++) {
            Rand[j * m + i].r =  RealMat[j * m + i];
            Rand[j * m + i].i =  0.0;
        }
    }

    /* Allocate work memory for QR.    */
    if (!(tau = (DoubleComplex *)
                    malloc (dim * sizeof (DoubleComplex)))) {
        printf ("randOrthMat: Malloc for *tau failed.\n");
        return 1;
    }
    if (!(work = (DoubleComplex *)
                    malloc (lwork * sizeof (DoubleComplex)))) {
        printf ("randOrthMat: Malloc for *work failed.\n");
        return 1;
    }

    /* QR factorize matrix Rand to get matrix of orthonormal columns.*/ 
    zgeqrf_ (&m, &n, Rand, &m, tau, work, &lwork, &info);
    zungqr_ (&m, &n, &dim, Rand, &m, tau, work, &lwork, &info);
    zlacpy_ ("", &m, &n, Rand, &m, A_m, &m);

    free (Rand);
    free (RealMat);
    free (tau);
    free (work);

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Generate a random vector. The real and imaginary part of each 
 *   element are in uniform (-1,1) distribution.  
 *
 * Input    
 *   x      Vector x
 *
 * Output     
 *   x      Random vector
 *
 ***************************************************************************/

void randComplexVec (DoubleComplexVec *x)
{
    zlarnv_ (&idist, &iseed[0] ,&(x->n), x->vec);
}

/***************************************************************************
 *
 * Abstract 
 *   Generate a random number.
 *
 * Input
 *   None
 *
 * Output
 *   None
 * 
 * Return
 *   The random number.
 *
 ***************************************************************************/

DoubleComplex randComplexNum ()
{
    DoubleComplex rand;
    int len = 1;
    
    zlarnv_ (&idist, &iseed[0] ,&len, &rand);

    return rand;
}