===================================================================
RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v
retrieving revision 1.133
retrieving revision 1.140
diff -u -p -r1.133 -r1.140
--- OpenXM_contrib2/asir2000/engine/nd.c	2006/06/05 08:11:10	1.133
+++ OpenXM_contrib2/asir2000/engine/nd.c	2006/06/12 21:51:33	1.140
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.132 2006/06/05 01:29:24 noro Exp $ */
+/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.139 2006/06/12 00:46:48 noro Exp $ */
 
 #include "nd.h"
 
@@ -1592,6 +1592,7 @@ NODE nd_gb(int m,int ishomo,int checkonly)
 	NDV nfv;
 	Q q,num,den;
 	union oNDC dn;
+	int diag_count = 0;
 
 	g = 0; d = 0;
 	for ( i = 0; i < nd_psn; i++ ) {
@@ -1604,6 +1605,7 @@ again:
 		l = nd_minp(d,&d);
 		if ( SG(l) != sugar ) {
 			if ( ishomo ) {
+				diag_count = 0;
 				stat = do_diagonalize(sugar,m);
 				if ( !stat ) {
 					NEXT(l) = d; d = l;
@@ -1639,6 +1641,15 @@ again:
 			}
 			nfv = ndtondv(m,nf); nd_free(nf);
 			nh = ndv_newps(m,nfv,0);
+			if ( ishomo && ++diag_count == diag_period ) {
+				diag_count = 0;
+				stat = do_diagonalize(sugar,m);
+				if ( !stat ) {
+					NEXT(l) = d; d = l;
+					d = nd_reconstruct(1,d);
+					goto again;
+				}
+			}
 			d = update_pairs(d,g,nh);
 			g = update_base(g,nh);
 			FREENDP(l);
@@ -2208,9 +2219,15 @@ void ndv_setup(int mod,int trace,NODE f,int dont_sort)
 	for ( nd_psn = 0, s = f; s; s = NEXT(s) ) if ( BDY(s) ) nd_psn++;
 	w = (NDV *)ALLOCA(nd_psn*sizeof(NDV));
 	for ( i = 0, s = f; s; s = NEXT(s) ) if ( BDY(s) ) w[i++] = BDY(s);
-	if ( !dont_sort )
-		qsort(w,nd_psn,sizeof(NDV),
-			(int (*)(const void *,const void *))ndv_compare);
+	if ( !dont_sort ) {
+		/* XXX heuristic */
+		if ( !nd_ord->id && (nd_ord->ord.simple<2) )
+			qsort(w,nd_psn,sizeof(NDV),
+				(int (*)(const void *,const void *))ndv_compare_rev);
+		else
+			qsort(w,nd_psn,sizeof(NDV),
+				(int (*)(const void *,const void *))ndv_compare);
+	}
 	nd_pslen = 2*nd_psn;
 	nd_ps = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
 	nd_ps_trace = (NDV *)MALLOC(nd_pslen*sizeof(NDV));
@@ -4168,9 +4185,26 @@ IndArray nm_ind_pair_to_vect_compress(int mod,UINT *s0
 	return r;
 }
 
+int compress_array(Q *svect,Q *cvect,int n)
+{
+	int i,j;
+
+	for ( i = j = 0; i < n; i++ )
+		if ( svect[i] ) cvect[j++] = svect[i];
+	return j;
+}
+
+void expand_array(Q *svect,Q *cvect,int n)
+{
+	int i,j;
+
+	for ( i = j = 0; j < n;  i++  )
+		if ( svect[i] ) svect[i] = cvect[j++];
+}
+
 int ndv_reduce_vect_q(Q *svect,int trace,int col,IndArray *imat,NM_ind_pair *rp0,int nred)
 {
-	int i,j,k,len,pos,prev;
+	int i,j,k,len,pos,prev,nz;
 	Q cs,mcs,c1,c2,cr,gcd,t;
 	IndArray ivect;
 	unsigned char *ivc;
@@ -4181,11 +4215,13 @@ int ndv_reduce_vect_q(Q *svect,int trace,int col,IndAr
 	NODE rp;
 	int maxrs;
 	double hmag;
-	struct oVECT v;
+	Q *cvect;
 
-	v.id = O_VECT; v.len = col; v.body = (pointer *)svect;
 	maxrs = 0;
-	hmag = p_mag((P)svect[0])*nd_scale;
+	for ( i = 0; i < col && !svect[i]; i++ );
+	if ( i == col ) return maxrs;
+	hmag = p_mag((P)svect[i])*nd_scale;
+	cvect = (Q *)ALLOCA(col*sizeof(Q));
 	for ( i = 0; i < nred; i++ ) {
 		ivect = imat[i];
 		k = ivect->head;
@@ -4224,11 +4260,19 @@ int ndv_reduce_vect_q(Q *svect,int trace,int col,IndAr
 					}
 					break;
 			}
+			for ( j = k+1; j < col && !svect[j]; j++ );
+			if ( j == col ) break;
+			if ( hmag && ((double)p_mag((P)svect[j]) > hmag) ) {
+				nz = compress_array(svect,cvect,col);
+				removecont_array(cvect,nz);
+				expand_array(svect,cvect,nz);
+				hmag = ((double)p_mag((P)svect[j]))*nd_scale;
+			}
 		}
-		if ( hmag && ((double)p_mag((P)svect[0]) > hmag) )
-			igcdv(&v,&t);
 	}
-	igcdv(&v,&t);
+	nz = compress_array(svect,cvect,col);
+	removecont_array(cvect,nz);
+	expand_array(svect,cvect,nz);
 	if ( DP_Print ) { 
 		fprintf(asir_out,"-"); fflush(asir_out);
 	}
@@ -4807,6 +4851,7 @@ NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,int trace,U
 	struct oEGT eg0,eg1,eg2,eg_f4,eg_f4_1,eg_f4_2;
 	int maxrs;
 	int *spsugar;
+	pointer *w;
 
 	spcol = col-nred;
 	get_eg(&eg0);
@@ -4840,12 +4885,19 @@ NODE nd_f4_red_q_main(ND_pairs sp0,int nsp,int trace,U
 	/* elimination (2nd step) */
 	colstat = (int *)ALLOCA(spcol*sizeof(int));
 	rank = nd_gauss_elim_q(spmat,spsugar,sprow,spcol,colstat);
-	r0 = 0;
+	w = (pointer *)ALLOCA(rank*sizeof(pointer));
 	for ( i = 0; i < rank; i++ ) {
-		NEXTNODE(r0,r); BDY(r) = 
-			(pointer)vect_to_ndv_q(spmat[i],spcol,col,rhead,s0vect);
-		SG((NDV)BDY(r)) = spsugar[i];
+		w[rank-i-1] = (pointer)vect_to_ndv_q(spmat[i],spcol,col,rhead,s0vect);
+		SG((NDV)w[rank-i-1]) = spsugar[i];
 /*		GC_free(spmat[i]); */
+	}
+#if 0
+	qsort(w,rank,sizeof(NDV),
+		(int (*)(const void *,const void *))ndv_compare);
+#endif
+	r0 = 0;
+	for ( i = 0; i < rank; i++ ) {
+		NEXTNODE(r0,r); BDY(r) = w[i];
 	}
 	if ( r0 ) NEXT(r) = 0;