Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
core/vgauss.c
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
113 lines (105 sloc)
3.8 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| /* Copyright 1990-2011, Jsoftware Inc. All rights reserved. */ | |
| /* License in license.txt. */ | |
| /* */ | |
| /* Verbs: Gaussian Elimination */ | |
| #include "j.h" | |
| F1(jtgausselm){A t;C*tv;I c,e,i,j,m,old,r,r1,*s;Q p,*u,*v,*x; | |
| F1RANK(2,jtgausselm,0); | |
| ASSERT(RAT&AT(w),EVNONCE); | |
| ASSERT(2==AR(w),EVRANK); | |
| s=AS(w); r=s[0]; c=s[1]; r1=MIN(r,c); v=QAV(w); m=c*bp(AT(w)); | |
| GA(t,LIT,m,1,0); tv=CAV(t); | |
| old=jt->tbase+jt->ttop; | |
| for(j=0;j<r1;++j){ | |
| e=-1; u=v+c*j+j; DO(r-j, if(XDIG(u->n)){e=i+j; break;} u+=c;); /* find pivot row */ | |
| ASSERT(0<=e,EVDOMAIN); | |
| x=v+c*j; | |
| if(j!=e){u=v+c*e; MC(tv,x,m); MC(x,u,m); MC(u,tv,m);} /* interchange rows e and j */ | |
| p=x[j]; DO(c, x[i]=qdiv(x[i],p);); | |
| for(i=0;i<r;++i){ | |
| if(i==j)continue; | |
| u=v+c*i; p=u[j]; /* pivot */ | |
| DO(c, u[i]=qminus(u[i],qtymes(p,x[i]));); | |
| } | |
| gc(w,old); | |
| } | |
| R w; | |
| } /* Gaussian elimination in place */ | |
| static F1(jtdetr){A t,z;C*tv;I c,e,g=1,i,j,k,m,old,r,*s;Q d,p,*u,*v,*x; | |
| RZ(w); | |
| s=AS(w); r=s[0]; c=s[1]; | |
| v=QAV(w); | |
| m=c*sizeof(Q); GA(t,LIT,m,1,0); tv=CAV(t); | |
| old=jt->tbase+jt->ttop; | |
| for(j=0;j<r;++j){ | |
| e=-1; u=v+c*j+j; DO(r-j, if(XDIG(u->n)){e=i+j; break;} u+=c;); /* find pivot row */ | |
| if(0>e)R cvt(RAT,zero); | |
| x=v+c*j; | |
| if(j!=e){u=v+c*e; MC(tv,x,m); MC(x,u,m); MC(u,tv,m); g=-g;} /* interchange rows e and j */ | |
| i=XDIG(x[j].n); if(i==XPINF||i==XNINF)R mark; | |
| for(i=j+1;i<r;++i){ | |
| u=v+c*i; | |
| if(XDIG(u[j].n)){p=qdiv(u[j],x[j]); for(k=j+1;k<r;++k)u[k]=qminus(u[k],qtymes(p,x[k]));} | |
| } | |
| gc(w,old); | |
| } | |
| d=0<g?*v:qminus(zeroQ,*v); u=v+1+c; DO(r-1, d=qtymes(d,*u); u+=1+c;); | |
| RE(0); | |
| GA(z,RAT,1,0,0); *QAV(z)=d; R z; | |
| } /* determinant on rational matrix; works in place */ | |
| static F1(jtdetd){D g,h,p,q,*u,*v,*x,*y,z=1.0;I c,d,e,i,j,k,r,*s; | |
| RZ(w); | |
| s=AS(w); r=s[0]; c=s[1]; v=DAV(w); | |
| NAN0; | |
| for(j=0;j<r;++j){ | |
| x=v+c*j; u=x+j; h=0.0; | |
| DO(r-j, k=i; DO(c-j, g=ABS(*u); if(h<g){h=g; d=j+k; e=j+i;} ++u;); u+=j;); /* find pivot, maximum abs element */ | |
| if(h==inf)R mark; | |
| if(0==h)R scf(0.0); | |
| if(j!=d){u=v+c*d+j; y=x+j; DO(c-j, q=*u; *u=*y; *y=q; ++u; ++y; ); z=-z;} /* interchange rows j and d */ | |
| if(j!=e){u=x+e; y=x+j; DO(r-j, q=*u; *u=*y; *y=q; u+=c; y+=c;); z=-z;} /* interchange cols j and e */ | |
| q=x[j]; z*=q; JBREAK0; | |
| for(i=j+1;i<r;++i){ | |
| u=v+c*i; | |
| if(u[j]){p=u[j]/q; for(k=j+1;k<r;++k)u[k]-=p*x[k];} | |
| }} | |
| NAN1; | |
| R scf(z); | |
| } /* determinant on real matrix; works in place */ | |
| #define ZABT(v) ((v).re*(v).re+(v).im*(v).im) | |
| static F1(jtdetz){A t;D g,h;I c,d,e,i,j,k,r,*s;Z p,q,*u,*v,*x,*y,z; | |
| RZ(w); | |
| z.re=1.0; z.im=0.0; | |
| s=AS(w); r=s[0]; c=s[1]; v=ZAV(w); | |
| NAN0; | |
| for(j=0;j<r;++j){ | |
| x=v+c*j; u=x+j; h=0.0; | |
| DO(r-j, k=i; DO(c-j, g=ZABT(*u); if(h<g){h=g; d=j+k; e=j+i;} ++u;); u+=j;); /* find pivot, maximum abs element */ | |
| if(h==inf)R mark; | |
| if(0==h)R scf(0.0); | |
| if(j!=d){u=v+c*d; DO(c-j, q=u[j+i]; u[j+i]=x[j+i]; x[j+i]=q;); z=zminus(zeroZ,z);} /* interchange rows j and d */ | |
| if(j!=e){u=v+c*j+e; y=v+c*j+j; DO(r-j, q=*u; *u=*y; *y=q; u+=c; y+=c;); z=zminus(zeroZ,z);} /* interchange cols j and e */ | |
| q=x[j]; z=ztymes(z,q); | |
| for(i=j+1;i<r;++i){ | |
| u=v+c*i; | |
| if(ZNZ(u[j])){p=zdiv(u[j],q); for(k=j+1;k<r;++k)u[k]=zminus(u[k],ztymes(p,x[k]));} | |
| }} | |
| NAN1; RE(0); | |
| GA(t,CMPX,1,0,0); *ZAV(t)=z; R t; | |
| } /* determinant on complex matrix; works in place */ | |
| F1(jtgaussdet){A z;I*s; | |
| RZ(w); | |
| ASSERT(2==AR(w),EVRANK); | |
| s=AS(w); | |
| ASSERT(s[0]==s[1],EVLENGTH); | |
| switch(AT(w)){ | |
| default: ASSERT(0,EVDOMAIN); | |
| case B01: | |
| case INT: R detd(cvt(FL,w)); | |
| case FL: z=detd(ca(w)); break; | |
| case CMPX: z=detz(ca(w)); break; | |
| case XNUM: z=detr(cvt(RAT,w)); break; | |
| case RAT: z=detr(ca(w)); | |
| } | |
| R z==mark?detxm(w,eval("-/ .*")):z; | |
| } /* determinant on square matrix */ |