Initial import

git-svn-id: http://moon:8086/svn/software/trunk/libsrc/_to_be_added@1 b431acfa-c32f-4a4a-93f1-934dc6c82436

mutil.c

@@ -0,0 +1,220 @@
+// ------------------------------------------------------------
+// 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;
+}