===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v
retrieving revision 1.46
retrieving revision 1.49
diff -u -p -r1.46 -r1.49
--- OpenXM_contrib2/asir2000/builtin/array.c	2005/02/08 18:06:05	1.46
+++ OpenXM_contrib2/asir2000/builtin/array.c	2005/12/21 23:18:15	1.49
@@ -45,7 +45,7 @@
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *
- * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.45 2005/01/23 14:03:47 noro Exp $
+ * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.48 2005/11/27 05:37:53 noro Exp $
 */
 #include "ca.h"
 #include "base.h"
@@ -64,7 +64,7 @@ extern int DP_Print; /* XXX */
 
 void Pnewvect(), Pnewmat(), Psepvect(), Psize(), Pdet(), Pleqm(), Pleqm1(), Pgeninvm();
 void Pinvmat();
-void Pnewbytearray();
+void Pnewbytearray(),Pmemoryplot_to_coord();
 
 void Pgeneric_gauss_elim();
 void Pgeneric_gauss_elim_mod();
@@ -107,6 +107,7 @@ struct ftab array_tab[] = {
 	{"matr",Pmat,-99999999},
 	{"matc",Pmatc,-99999999},
 	{"newbytearray",Pnewbytearray,-2},
+	{"memoryplot_to_coord",Pmemoryplot_to_coord,1},
 	{"sepmat_destructive",Psepmat_destructive,2},
 	{"sepvect",Psepvect,2},
 	{"qsort",Pqsort,-2},
@@ -477,6 +478,37 @@ void Pnewbytearray(NODE arg,BYTEARRAY *rp)
 	*rp = array;
 }
 
+#define MEMORY_GETPOINT(a,len,x,y) (((a)[(len)*(y)+((x)>>3)])&(1<<((x)&7)))
+
+void Pmemoryplot_to_coord(NODE arg,LIST *rp)
+{
+	int len,blen,y,i,j;
+	char *a;
+	NODE r0,r,n;
+	LIST l;
+	BYTEARRAY ba;
+	Q iq,jq;
+
+	asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");
+	arg = BDY((LIST)ARG0(arg));
+	len = QTOS((Q)ARG0(arg));
+	blen = (len+7)/8;
+	y = QTOS((Q)ARG1(arg));
+	ba = (BYTEARRAY)ARG2(arg); a = ba->body;
+	r0 = 0;
+	for ( j = 0; j < y; j++ )
+		for ( i = 0; i < len; i++ )
+			if ( MEMORY_GETPOINT(a,blen,i,j) ) {
+				NEXTNODE(r0,r);
+				STOQ(i,iq); STOQ(j,jq);
+				n = mknode(2,iq,jq);
+				MKLIST(l,n);
+				BDY(r) = l;
+			}
+	if ( r0 ) NEXT(r) = 0;
+	MKLIST(*rp,r0);
+}
+
 void Pnewmat(NODE arg,MAT *rp)
 {
 	int row,col;
@@ -842,17 +874,34 @@ void Pinvmat(NODE arg,LIST *rp)
 
 void Pgeneric_gauss_elim(NODE arg,LIST *rp)
 {
-	NODE n0;
+	NODE n0,opt,p;
 	MAT m,nm;
 	int *ri,*ci;
 	VECT rind,cind;
 	Q dn,q;
 	int i,j,k,l,row,col,t,rank;
+	int is_hensel = 0;
+	char *key;
+	Obj value;
 
+	if ( current_option ) {
+		for ( opt = current_option; opt; opt = NEXT(opt) ) {
+			p = BDY((LIST)BDY(opt));
+			key = BDY((STRING)BDY(p));
+			value = (Obj)BDY(NEXT(p));
+			if ( !strcmp(key,"hensel") && value ) {
+				is_hensel = value ? 1 : 0;
+				break;
+			}
+		}
+	}
 	asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");
 	m = (MAT)ARG0(arg);
 	row = m->row; col = m->col;
-	rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
+	if ( is_hensel )
+		rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);
+	else
+		rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
 	t = col-rank;
 	MKVECT(rind,rank);
 	MKVECT(cind,t);
@@ -876,6 +925,8 @@ void Pgeneric_gauss_elim(NODE arg,LIST *rp)
 		B : a rank(A) x col-rank(A) matrix
 		R : a vector of length rank(A)
 		C : a vector of length col-rank(A)
+		RN : a vector of length rank(A) indicating useful rows
+
 		B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
 */
 
@@ -883,11 +934,11 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 {
 	NODE n0;
 	MAT m,mat;
-	VECT rind,cind;
+	VECT rind,cind,rnum;
 	Q **tmat;
-	int **wmat;
-	Q *rib,*cib;
-	int *colstat;
+	int **wmat,**row0;
+	Q *rib,*cib,*rnb;
+	int *colstat,*p;
 	Q q;
 	int md,i,j,k,l,row,col,t,rank;
 
@@ -896,6 +947,10 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 	m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
 	row = m->row; col = m->col; tmat = (Q **)m->body;
 	wmat = (int **)almat(row,col);
+
+	row0 = (int **)ALLOCA(row*sizeof(int *));
+	for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
+
 	colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
 	for ( i = 0; i < row; i++ )
 		for ( j = 0; j < col; j++ )
@@ -908,6 +963,13 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 				wmat[i][j] = 0;
 	rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
 
+	MKVECT(rnum,rank);
+	rnb = (Q *)rnum->body;
+	for ( i = 0; i < rank; i++ )
+		for ( j = 0, p = wmat[i]; j < row; j++ )
+			if ( p == row0[j] )
+				STOQ(j,rnb[i]);
+
 	MKMAT(mat,rank,col-rank);
 	tmat = (Q **)mat->body;
 	for ( i = 0; i < rank; i++ )
@@ -925,7 +987,7 @@ void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 		} else {
 			STOQ(j,cib[l]); l++;
 		}
-	n0 = mknode(3,mat,rind,cind);
+	n0 = mknode(4,mat,rind,cind,rnum);
 	MKLIST(*rp,n0);
 }
 
@@ -1200,7 +1262,13 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 				} else
 					wi[j] = 0;
 
+		if ( DP_Print ) {
+			fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
+		}
 		rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
+		if ( DP_Print ) {
+			fprintf(asir_out,"done.\n"); fflush(asir_out);
+		}
 		a = (Q **)almat_pointer(rank,rank); /* lhs mat */
 		MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
 		for ( j = li = ri = 0; j < col; j++ )
@@ -1237,7 +1305,7 @@ int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn
 			wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
 			for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
 			for ( q = ONE, count = 0; ; ) {
-				if ( DP_Print > 3 )
+				if ( DP_Print )
 					fprintf(stderr,"o");
 				/* wc = -b mod md */
 				get_eg(&tmp0);