=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/builtin/array.c,v retrieving revision 1.47 retrieving revision 1.49 diff -u -p -r1.47 -r1.49 --- OpenXM_contrib2/asir2000/builtin/array.c 2005/11/27 00:07:05 1.47 +++ 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.46 2005/02/08 18:06:05 saito 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); @@ -1213,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++ ) @@ -1250,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);