/****************************************************************************
 *
 * MatCompute.c
 *
 * Abstract 
 *   Implement computing routines which are required by the block Lanczos
 *   algorithm, but are not provided by LAPACK or BLAS explicitly.
 *
 * Student:  Guohong Liu
 *
 * ID:       0385117
 *
 ****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include "computing.h"

/***************************************************************************
 *
 * Abstract 
 *   Compute the element-by-element product of two matrices,
 *      R = alpha * A .* B + C
 *   The implementation of it is referred to BLAS zaxpy_ ()
 * 
 * Input    
 *   alpha  Real constant.
 *   A      Matrix. Size = m * n
 *   B      Matrix. Size = m * n
 *   C      Matrix. Size = m * n
 *           
 * Output   
 *   R      Matrix. Size = m * n.
 *
 * Return  
 *   0      Successful exit
 *   1      Exception occurs
 *
 ***************************************************************************/

int memult (DoubleReal alpha, DoubleComplexMat *A, DoubleComplexMat *B, 
            DoubleComplexMat *C, DoubleComplexMat *R)
{
    int i, n;
    DoubleComplex *x, *y, *z, *r, tmp;

    x = A->mat;
    y = B->mat;
    if (C != NULL) {
        z = C->mat;
    }
    r = R->mat;

    n = (A->m) * (A->n);
    if (n <= 0) {
	    return 1;
    }

    if (C != NULL) {
        for (i = 0; i < n; i++) {
	        tmp.r = alpha * x[i].r;
            tmp.i = alpha * x[i].i;
	        r[i].r = y[i].r * tmp.r - y[i].i * tmp.i + z[i].r;
            r[i].i = y[i].r * tmp.i + y[i].i * tmp.r + z[i].i;
        }
    } 
    else {
        for (i = 0; i < n; i++) {
	        tmp.r = alpha * x[i].r;
            tmp.i = alpha * x[i].i;
	        r[i].r = y[i].r * tmp.r - y[i].i * tmp.i;
            r[i].i = y[i].r * tmp.i + y[i].i * tmp.r;
        }
    }

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute the element-by-element product of two vectors,
 *      r = x .* y + z
 *      r = x .* conjg(y) + z
 *   where x is a complex vector. The implementation of it is referred 
 *   to BLAS zaxpy_ ().
 *
 * Input    
 *    x       Vector of order n
 *    y       Vector of order n.
 *    transy  "N", op(y) = y
 *            "C", op(y) = conjg(y)
 *    z       Vector.
 *
 * Output   
 *    r       Result vector
 *
 * Return  
 *    0       Successful exit
 *    1       Exception occurs
 *
 ***************************************************************************/

int vzemult (DoubleComplexVec *x, DoubleComplexVec *y, char *transy, 
             DoubleComplexVec *z, DoubleComplexVec *r)
{
    int i, n;
    DoubleComplex *x_m, *y_m, *z_m, *r_m;

    x_m = x->vec;
    y_m = y->vec;
    if (z != NULL) {
        z_m = z->vec;
    }
    r_m = r->vec;

    n = x->n;
    if (n <= 0) {
	    return 1;
    }

    if (transy[0] == 'N') {
        if (z != NULL) {
            for (i = 0; i < n; i++) {
	            r_m[i].r = y_m[i].r * x_m[i].r 
                           - y_m[i].i * x_m[i].i + z_m[i].r;
                r_m[i].i = y_m[i].r * x_m[i].i 
                           + y_m[i].i * x_m[i].r + z_m[i].i;
            }
        } 
        else {
            for (i = 0; i < n; i++) {
	            r_m[i].r = y_m[i].r * x_m[i].r - y_m[i].i * x_m[i].i;
                r_m[i].i = y_m[i].r * x_m[i].i + y_m[i].i * x_m[i].r;
            }
        }
    }
    else if (transy[0] == 'C') {
        if (z != NULL) {
            for (i = 0; i < n; i++) {
	            r_m[i].r = y_m[i].r * x_m[i].r 
                           + y_m[i].i * x_m[i].i + z_m[i].r;
                r_m[i].i = y_m[i].r * x_m[i].i 
                           - y_m[i].i * x_m[i].r + z_m[i].i;
            }
        }
        else {
            for (i = 0; i < n; i++) {
	            r_m[i].r = y_m[i].r * x_m[i].r + y_m[i].i * x_m[i].i;
                r_m[i].i = y_m[i].r * x_m[i].i - y_m[i].i * x_m[i].r;
            }
        }
    }

    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute the element-by-element product of two vectors,
 *      r = alpha * x .* y + z
 *      r = alpha * x .* conjg(y) + z
 *   where alpha is real number, and x is a real vector. The 
 *   implementation of it is referred to BLAS zaxpy_ ().
 *
 * Input    
 *    alpha   Real number.
 *    x       Vector of order n
 *    y       Vector of order n.
 *    transy  "N", op(y) = y
 *            "C", op(y) = conjg(y)
 *    z       Vector.
 *
 * Output   
 *    r       Result vector
 *
 * Return  
 *    0       Successful exit
 *    1       Exception occurs
 *
 ***************************************************************************/

int vdemult (DoubleReal alpha, DoubleRealVec *x, 
             DoubleComplexVec *y, char *transy, 
             DoubleComplexVec *z, DoubleComplexVec *r)
{
    int i, n;
    DoubleComplex *y_m, *z_m, *r_m;
    DoubleReal *x_m, tmp;

    x_m = x->vec;
    y_m = y->vec;
    if (z != NULL) {
        z_m = z->vec;
    }
    r_m = r->vec;

    n = x->n;
    if (n <= 0) {
	    return 1;
    }

    if (transy[0] == 'N') {
        if (z != NULL) {
            for (i = 0; i < n; i++) {
                tmp = alpha * x_m[i];
	            r_m[i].r = y_m[i].r * tmp + z_m[i].r;
                r_m[i].i = y_m[i].i * tmp + z_m[i].i;
            }
        } 
        else {
            for (i = 0; i < n; i++) {
                tmp = alpha * x_m[i];
	            r_m[i].r = y_m[i].r * tmp;
                r_m[i].i = y_m[i].i * tmp;
            }
        }
    }
    else if (transy[0] == 'C') {
        if (z != NULL) {
            for (i = 0; i < n; i++) {
                tmp = alpha * x_m[i];
	            r_m[i].r = y_m[i].r * tmp + z_m[i].r;
                r_m[i].i = -y_m[i].i * tmp + z_m[i].i;
            }
        }
        else {
            for (i = 0; i < n; i++) {
                tmp = alpha * x_m[i];
	            r_m[i].r = y_m[i].r * tmp;
                r_m[i].i = -y_m[i].i * tmp;
            }
        }
    }

    return 0;
}

/***************************************************************************
 *
 * Abstract
 *   Compute the multiplication and addition of three complex numbers,
 *      r =  a * b + c
 *      r = -a * b + c
 *      r =  a * conjg(b) + c
 *      r = -a * conjg(b) + c
 *   wher a, b, and c are complex numbers.
 *
 * Input    
 *   signa   "+" or "-"
 *   a       Complex number.
 *   b       Complex number.
 *   transb  "N", op(b) = b 
 *           "C", op(b) = conjg (b)   conjugate of b.
 *   c       Complex number.
 *
 * Output     
 *   r       Complex number.
 *
 ***************************************************************************/

int zzmult (char *signa, DoubleComplex *a, DoubleComplex *b, char *transb, 
            DoubleComplex *c, DoubleComplex *r)
{
    DoubleComplex result, b1;

    if (transb[0] == 'N') {
        b1.r = b->r;
        b1.i = b->i;
    }
    else if (transb[0] == 'C') {
        b1.r = b->r;
        b1.i = -b->i;
    }
    else {
        printf ("zmult: Wrong operation of b.\n");
        return 1;
    }

    if (signa[0] == '+') {
        result.r = a->r * b1.r - a->i * b1.i; 
        result.i = a->r * b1.i + a->i * b1.r;
    }
    else if (signa[0] == '-') {
        result.r = (-a->r) * b1.r + a->i * b1.i; 
        result.i = (-a->r) * b1.i - a->i * b1.r;    
    }
    else {
        printf ("zmult: Wrong sign of a.\n");
        return 1;
    }
    if (c != NULL) {
        r->r = result.r + c->r;
        r->i = result.i + c->i;
    }
    else {
        r->r = result.r;
        r->i = result.i;
    }
    
    return 0;
}

/***************************************************************************
 *
 * Abstract
 *   Compute the multiplication and addition of three complex numbers,
 *      r =  a * b + c
 *      r = -a * b + c
 *      r =  a * conjg(b) + c
 *      r = -a * conjg(b) + c
 *   where a is a real, b and c are complex numbers.
 *
 * Input    
 *   signa   "+" or "-"
 *   a       Real number.
 *   b       Complex number.
 *   transb  "N", op(b) = b 
 *           "C", op(b) = conjg (b)   conjugate of b.
 *   c       Complex number.
 *
 * Output     
 *   r       Complex number.
 *
 ***************************************************************************/

int dzmult (char *signa, DoubleReal a, DoubleComplex *b, char *transb, 
            DoubleComplex *c, DoubleComplex *r)
{
    DoubleComplex result, b1;

    if (transb[0] == 'N') {
        b1.r = b->r;
        b1.i = b->i;
    }
    else if (transb[0] == 'C') {
        b1.r = b->r;
        b1.i = -b->i;
    }
    else {
        printf ("zmult: Wrong operation of b.\n");
        return 1;
    }

    if (signa[0] == '+') {
        result.r = a * b1.r; 
        result.i = a * b1.i;
    }
    else if (signa[0] == '-') {
        result.r = -a * b1.r; 
        result.i = -a * b1.i;    
    }
    else {
        printf ("zmult: Wrong sign of a.\n");
        return 1;
    }
    if (c != NULL) {
        r->r = result.r + c->r;
        r->i = result.i + c->i;
    }
    else {
        r->r = result.r;
        r->i = result.i;
    }
    
    return 0;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute r = a / b, where a is a complex number and b is a real number.
 *
 * Input    
 *   a      Vector a
 *   b      Vector b
 *
 * Output     
 *   r      Result vector
 *
 ***************************************************************************/

void zdiv (DoubleComplex *a, DoubleReal b, DoubleComplex *r)
{
    r->r = a->r / b;
    r->i = a->i / b;
}

/***************************************************************************
 *
 * Abstract 
 *   Compute the conjugate of a complex number.
 *
 * Input    
 *   z      complex number
 *
 * Output     
 *   None
 *
 * Return
 *   the conjugate of z
 *
 ***************************************************************************/

DoubleComplex conjg (DoubleComplex z)
{
    DoubleComplex x;

    x.r = z.r;
    x.i = -z.i;

    return x;
}