Now wrangled with ability to use submatricies.

3 files changed, 55 insertions, 43 deletions

diff --git a/redist/dclapack.h b/redist/dclapack.h
index 76ce545..1dc7b77 100644
--- a/redist/dclapack.h
+++ b/redist/dclapack.h
@@ -10,16 +10,14 @@
 #define _ABS(a)  ( (a)<=0 ? 0-(a) : (a) )
 
 #ifdef DYNAMIC_INDEX
-	#define _(A,O,P) A[O*n+P]
-	#define _I(A,O,P,n) A[O*n+P]
+	#define _(A,O,P) A[O*A##c+P]
 #else
 	#define _(A,O,P) A[O][P]
-	#define _I(A,O,P,n) A[O][P]
 #endif
 
 
 /*
- * Prints a matrix A (n by m)
+ * Prints a matrix A (n by m)  (with witdth Ac)
  */
 #define PRINT(A,n,m) { \
     int i,j; \
@@ -34,7 +32,7 @@
 } 
 
 /*
- * Returns the identity matrix
+ * Returns the identity matrix (with size n x n, but width Ic)
  */
 #define IDENTITY(I,n) { \
     int i,j; \
@@ -148,11 +146,11 @@
 
 #ifdef DYNAMIC_INDEX
 	#define INV_SETUP(ORDER) \
-    FLOAT Ipiv[ORDER*ORDER]; \
-    FLOAT L[ORDER*ORDER]; \
-    FLOAT U[ORDER*ORDER]; \
-    FLOAT I[ORDER*ORDER]; \
-    FLOAT C[ORDER*ORDER];
+    FLOAT Ipiv[ORDER*ORDER]; const int Ipivc = ORDER; \
+    FLOAT L[ORDER*ORDER]; const int Lc = ORDER; \
+    FLOAT U[ORDER*ORDER]; const int Uc = ORDER; \
+    FLOAT I[ORDER*ORDER]; const int Ic = ORDER; \
+    FLOAT C[ORDER*ORDER]; const int Cc = ORDER;
 #else
 	#define INV_SETUP(ORDER) \
     FLOAT Ipiv[ORDER][ORDER]; \
@@ -194,7 +192,7 @@ PRINT(Ainv,n,n); \
         for (j=0; j<p; j++) { \
             _(R,i,j) = 0.0f; \
             for (k=0; k<m; k++) { \
-                _(R,i,j) += _(A,i,k) * _I(B,k,j,m); \
+                _(R,i,j) += _(A,i,k) * _(B,k,j); \
             } \
         } \
     } \
@@ -213,7 +211,7 @@ PRINT(Ainv,n,n); \
         for (j=0; j<p; j++) { \
             _(R,i,j) = _(C,i,j); \
             for (k=0; k<m; k++) { \
-                _(R,i,j) += _(A,i,k) * _I(B,k,j,m); \
+                _(R,i,j) += _(A,i,k) * _(B,k,j); \
             } \
         } \
     } \
@@ -233,7 +231,7 @@ PRINT(Ainv,n,n); \
         for (j=0; j<p; j++) { \
             sum = 0.0f; \
             for (k=0; k<m; k++) { \
-                sum += _(A,i,k) * _I(B,k,j,m); \
+                sum += _(A,i,k) * _(B,k,j); \
             } \
             _(R,i,j) = alpha * sum + beta * _(C,i,j); \
         } \
diff --git a/redist/dclhelpers.c b/redist/dclhelpers.c
index 23c3ca5..9fa9266 100644
--- a/redist/dclhelpers.c
+++ b/redist/dclhelpers.c
@@ -3,59 +3,59 @@
 #define DYNAMIC_INDEX
 #include <stdio.h>
 #include "dclapack.h"
+#include <alloca.h>
 
-
-void dclPrint( const DCL_FLOAT * A, int n, int m )
+void dclPrint( const DCL_FLOAT * A, int Ac, int n, int m )
 {
 	PRINT( A, n, m );
 }
 
-void dclIdentity( DCL_FLOAT * A, int n )
+void dclIdentity( DCL_FLOAT * I, int Ic, int n )
 {
-	IDENTITY( A, n );
+	IDENTITY( I, n );
 }
 
-void dclTransp( DCL_FLOAT * R, const DCL_FLOAT * A, int n, int m )
+void dclTransp( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int n, int m )
 {
 	TRANSP(R,A,n,m);
 }
 
-void dclLU( DCL_FLOAT * L, DCL_FLOAT * U, const DCL_FLOAT * A, int * Piv, int n )
+void dclLU( DCL_FLOAT * L, int Lc, DCL_FLOAT * U, int Uc, const DCL_FLOAT * A, int Ac, int * Piv, int n )
 {
 	LU(L,U,A,Piv,n);
 }
 
-void dclPivot( DCL_FLOAT * R, const DCL_FLOAT * A, int * Piv, int n, int m )
+void dclPivot( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int * Piv, int n, int m )
 {
 	PIVOT(R,A,Piv,n,m);
 }
 
-void dclLSub( DCL_FLOAT * X, const DCL_FLOAT * L, const DCL_FLOAT * B, int n, int m )
+void dclLSub( DCL_FLOAT * X, int Xc, const DCL_FLOAT * L, int Lc, const DCL_FLOAT * B, int Bc, int n, int m )
 {
 	L_SUB(X,L,B,n,m);
 }
 
-void dclUSub( DCL_FLOAT * X, const DCL_FLOAT * U, const DCL_FLOAT * B, int n, int m )
+void dclUSub( DCL_FLOAT * X, int Xc, const DCL_FLOAT * U, int Uc, const DCL_FLOAT * B, int Bc, int n, int m )
 {
 	U_SUB(X,U,B,n,m);
 }
 
-void dclInv( DCL_FLOAT * Ainv, const DCL_FLOAT * A, int n )
+void dclInv( DCL_FLOAT * Ainv, int Ainvc, const DCL_FLOAT * A, int Ac, int n )
 {
 	INV(Ainv,A,n,n);
 }
 
-void dclMul( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, int n, int m, int p )
+void dclMul( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, int n, int m, int p )
 {
 	MUL(R,A,B,n,m,p);
 }
 
-void dclMulAdd( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, int n, int m, int p )
+void dclMulAdd( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, const DCL_FLOAT * C, int Cc, int n, int m, int p )
 {
 	MULADD(R,A,B,C,n,m,p);
 }
 
-void dclGMulAdd( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p )
+void dclGMulAdd( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, const DCL_FLOAT * C, int Cc, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p )
 {
 	GMULADD(R,A,B,C,alpha,beta,n,m,p);
 }
@@ -77,12 +77,26 @@ void dcldgemm(
 	int ldc //must be n
 	 )
 {
-	DCL_FLOAT * ta;
-	DCL_FLOAT * tb;
+	const DCL_FLOAT * ta;
+	const DCL_FLOAT * tb;
 	if( transA )
 	{
 		ta = alloca( sizeof( DCL_FLOAT ) * n * m );
-		
+		//TRANSP(	ta, A, n, m );
+	}
+	else
+		ta = A;
+
+	if( transB )
+	{
+		tb = alloca( sizeof( DCL_FLOAT ) * n * m );
+		//TRANSP(	tb, B, n, m );
 	}
+	else
+		tb = B;
+
+
+	//XXX Tricky: In here, A is mxn
+
 }
 
diff --git a/redist/dclhelpers.h b/redist/dclhelpers.h
index fd4e02d..5010e02 100644
--- a/redist/dclhelpers.h
+++ b/redist/dclhelpers.h
@@ -6,56 +6,56 @@
 //XXX XXX XXX WARNING XXX XXX XXX The argument order may be changing!!!
 
 /* Prints matrix A of size[n][m] */
-void dclPrint( const DCL_FLOAT * A, int n, int m );
+void dclPrint( const DCL_FLOAT * A, int Ac, int n, int m );
 
 /* Returns the identity matrix */
-void dclIdentity( DCL_FLOAT * A, int n );
+void dclIdentity( DCL_FLOAT * I, int Ic, int n );
 
 /* R = Transpose(A)
      A is (n by m)
      R is (m by n) */
-void dclTransp( DCL_FLOAT * R, const DCL_FLOAT * A, int n, int m );
+void dclTransp( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int n, int m );
 
 /* Calculate L,U of a matrix A with pivot table; the pivot table is output. */
-void dclLU( DCL_FLOAT * L, DCL_FLOAT * U, const DCL_FLOAT * A, int * Piv, int n );
+void dclLU( DCL_FLOAT * L, int Lc, DCL_FLOAT * U, int Uc, const DCL_FLOAT * A, int Ac, int * Piv, int n );
 
 /* Pivots a matrix to a different matrix
     R = Pivot(A) given table 'Piv'
     A and R are (n by m) */
-void dclPivot( DCL_FLOAT * R, const DCL_FLOAT * A, int * Piv, int n, int m );
+void dclPivot( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, int * Piv, int n, int m );
 
 /* Solve LX=B for matrix X and B
     L is n by n (lower triangular)
     B is n by m */
-void dclLSub( DCL_FLOAT * X, const DCL_FLOAT * L, const DCL_FLOAT * B, int n, int m );
+void dclLSub( DCL_FLOAT * X, int Xc, const DCL_FLOAT * L, int Lc, const DCL_FLOAT * B, int Bc, int n, int m );
 
 /* Solve UX=B for matrix X and B
     U is n by n (upper triangular)
     B is n by m */
-void dclUSub( DCL_FLOAT * X, const DCL_FLOAT * U, const DCL_FLOAT * B, int n, int m );
+void dclUSub( DCL_FLOAT * X, int Xc, const DCL_FLOAT * U, int Uc, const DCL_FLOAT * B, int Bc, int n, int m );
 
 /* Inverts a matrix X (n by n) using the method of LU decomposition */
-void dclInv( DCL_FLOAT * Ainv, const DCL_FLOAT * A, int n );
+void dclInv( DCL_FLOAT * Ainv, int Ainvc, const DCL_FLOAT * A, int Ac, int n );
 
 /* Matrix Multiply R = A * B
     A (n by m)
     B (m by p)
     R (n by p) */
-void dclMul( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, int n, int m, int p );
+void dclMul( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, int n, int m, int p );
 
 /* Matrix Multiply R = A * B + C
     A (n by m)
     B (m by p)
     C (n by p)
     R (n by p) */
-void dclMulAdd( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, int n, int m, int p );
+void dclMulAdd( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, const DCL_FLOAT * C, int Cc, int n, int m, int p );
 
 /* Matrix Multiply R = alpha * A * B + beta * C
     A (n by m)
     B (m by p)
     C (n by p)
     R (n by p) */
-void dclGMulAdd( DCL_FLOAT * R, const DCL_FLOAT * A, const DCL_FLOAT * B, const DCL_FLOAT * C, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p );
+void dclGMulAdd( DCL_FLOAT * R, int Rc, const DCL_FLOAT * A, int Ac, const DCL_FLOAT * B, int Bc, const DCL_FLOAT * C, int Cc, DCL_FLOAT alpha, DCL_FLOAT beta, int n, int m, int p );
 
 
 /********************************
@@ -71,12 +71,12 @@ void dcldgemm(
 	int k,
 	DCL_FLOAT alpha,
 	const DCL_FLOAT* A,
-	int lda, //must be n
+	int Ac,
 	const DCL_FLOAT* B,
-	int ldb, //must be m
+	int Bc,
 	DCL_FLOAT beta,
 	const DCL_FLOAT * C,
-	int ldc //must be n
+	int Cc
 	 );