// ------------------------------------------------------------
// mutil.c
// Matrix utilities
// 5.3.2005, J.Ahrensfeld
//
// ------------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "lapack_cwrap.h"

// ------------------------------------------------------------
double lp_work_buf[1024] = {0};

// ------------------------------------------------------------
void Eye(double *pI, double value, int N)
{
	int len;

	memset(pI, 0, N*N*sizeof(double));

	len = 0;

	while(len < N*N)
	{
		pI[len] = value;
		len += (N+1);
	}

}

// ------------------------------------------------------------
void RandMat(double *pA, int N)
{
	int i;

	for(i=0; i < N; i++)
		pA[i] = 2.0*(0.5 - (double)rand()/RAND_MAX);
}

// ------------------------------------------------------------
void IncMat(double *pA, int nRow, int nCol, char order)
{
	int i, j, k;

	if (order == 'R')
	{
		k = 1;
		for(i=0; i < nCol; i++)
			for(j=0; j < nRow; j++)
				pA[i*nRow+j] = (double)k++;
	}
	else
	{
		k = 1;
		for(i=0; i < nRow; i++)
			for(j=0; j < nCol; j++)
				pA[j*nRow+i] = (double)k++;
	}
}

// ------------------------------------------------------------
void CopyTrans(double *pA, double *pAT, int nRow, int nCol)
{
	int i, j;

	for(i=0; i < nRow; i++)
		for(j=0; j < nCol; j++)
			pAT[j+i*nCol] = pA[j*nRow+i];

}

// ------------------------------------------------------------
void GetDiag(double *pA, double *pY, int N)
{
	int i;
	for (i=0; i < N; i++)
		pY[i] = pA[i*N+i];

}

// ------------------------------------------------------------
void MakeDiag(double *pD, double *pX, int N)
{
	int len;

	memset(pD, 0, N*N*sizeof(double));

	len = 0;
	while(len < N*N)
	{
		pD[len] = *(pX++);
		len += (N+1);
	}

}

// ------------------------------------------------------------
void PrintMatrix(char *pStr, double *A, int num_row, int num_col)
{
	int i, j;
	double v;

	printf("%s\n",pStr);
    for (i=0; i<num_row; i++)
	{
	    for (j=0; j<num_col; j++)
	    {
	    	v = A[j*num_row+i];
			printf ("%5.5g ", v);
		}

		printf("\n");
	}

}

// ------------------------------------------------------------
void PrintMatrixSym(char *pStr, double *A, int num_row, int num_col, char uplo)
{
	int i, j;
	double v;

	printf("%s\n",pStr);
    for (i=0; i<num_row; i++)
	{
	    for (j=0; j<num_col; j++)
	    {
	    	v = A[j*num_row+i];
			printf ("%5.5g ", v);
		}

		printf("\n");
	}

}

// ------------------------------------------------------------
double smin(double v1, double v2)
{
	if (v1 < v2)
		return v1;

	return v2;
}

// ------------------------------------------------------------
double smax(double v1, double v2)
{
	if (v1 > v2)
		return v1;

	return v2;
}

// ------------------------------------------------------------
void vsum(double *pX, double *pY, int lda, int N)
{
	int i, j;

	for (j=0; j < lda; j++)
		pY[j] = *(pX++);

	for(i=1; i < N; i++)
		for (j=0; j < lda; j++)
			pY[j] += *(pX++);

}

// ------------------------------------------------------------
void GetPrincipleEigenVector(double *pA, double *pV, int lda)
{
	double max_x;
	int max_i, i;

	max_x = pV[0];
	max_i = 0;

	for (i=1; i < lda; i++)
	{
		if(pV[i] > max_x)
		{
			max_x = pV[0];
			max_i = i;
		}
	}

	for (i=0; i < lda; i++)
		pV[i] = pA[max_i*lda + i];
}

// ------------------------------------------------------------
double dnrm(int N, int M, double *A, char mode)
{
	int i;
	double maxval, val;

	maxval = 0;
	if (mode == '1')
	{
		for (i=0; i < N*M; i += N)
		{
			val = dasum(N, &A[i], 1);
			if (val > maxval)
				maxval = val;
		}
	}
	else
	{
		for (i=0; i < N; i++)
		{
			val = dasum(M, &A[i], N);
			if (val > maxval)
				maxval = val;
		}
	}

	return maxval;
}