Mercurial > hg > jgplsrc

annotate vd.c @ 0:e0bbaa717f41 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
lol J
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Mon, 25 Nov 2013 11:56:30 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
1 /* Copyright 1990-2011, Jsoftware Inc. All rights reserved. */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
2 /* License in license.txt. */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
3 /* */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
4 /* Verbs: Domino */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
5
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
6 #include "j.h"
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
7
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
8
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
9 static F1(jtnorm){R sqroot(pdt(w,conjug(w)));}
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
10
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
11 F1(jtrinv){PROLOG;A ai,bx,di,z;I m,n,r,*s;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
12 F1RANK(2,jtrinv,0);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
13 r=AR(w); s=AS(w); n=2>r?1:s[1]; m=(1+n)/2;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
14 ASSERT(!r||n==s[0],EVLENGTH);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
15 if(1>=n)R recip(w);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
16 ai=rinv(take(v2(m,m),w));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
17 di=rinv(drop(v2(m,m),w));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
18 bx=negate(pdt(ai,pdt(take(v2(m,m-n),w),di)));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
19 z=over(stitch(ai,bx),take(v2(n-m,-n),di));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
20 EPILOG(z);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
21 } /* R.K.W. Hui, Uses of { and }, APL87, p. 56 */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
22
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
23 static F1(jtqrr){PROLOG;A a1,q,q0,q1,r,r0,r1,t,*tv,t0,t1,y,z;I m,n,p,*s;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
24 RZ(w);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
25 if(2>AR(w)){p=AN(w); n=m=1;}else{s=AS(w); p=s[0]; n=s[1]; m=(1+n)/2;}
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
26 if(1>=n){
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
27 t=norm(ravel(w));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
28 ASSERT(!AN(w)||!equ(t,zero),EVDOMAIN);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
29 RZ(q=divide(w,t));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
30 R link(2>AR(q)?table(q):q,reshape(v2(n,n),p?t:one));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
31 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
32 RZ(t0=qrr(take(v2(p,m),w)));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
33 tv=AAV(t0); q0=*tv++; r0=*tv;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
34 RZ(a1=drop(v2(0L,m),w));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
35 RZ(y=pdt(conjug(cant1(q0)),a1));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
36 RZ(t1=qrr(minus(a1,pdt(q0,y))));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
37 tv=AAV(t1); q1=*tv++; r1=*tv;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
38 RZ(q=stitch(q0,q1));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
39 RZ(r=over(stitch(r0,y),take(v2(n-m,-n),r1)));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
40 z=link(q,r); EPILOG(z);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
41 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
42
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
43 F1(jtqr){A r,z;D c=inf,d=0,x;I n1,n,*s,wr;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
44 F1RANK(2,jtqr,0);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
45 ASSERT(DENSE&AT(w),EVNONCE);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
46 ASSERT(AT(w)&B01+INT+FL+CMPX,EVDOMAIN);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
47 wr=AR(w); s=AS(w);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
48 ASSERT(2>wr||s[0]>=s[1],EVLENGTH);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
49 RZ(z=qrr(w)); r=*(1+AAV(z)); n=*AS(r); n1=1+n;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
50 if(FL&AT(r)){D*v=DAV(r); DO(n, x= ABS(*v); if(x<c)c=x; if(x>d)d=x; v+=n1;);}
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
51 else {Z*v=ZAV(r); DO(n, x=zmag(*v); if(x<c)c=x; if(x>d)d=x; v+=n1;);}
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
52 ASSERT(!n||c>d*jt->fuzz,EVDOMAIN);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
53 R z;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
54 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
55
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
56 static F2(jticor){D d,*v;I n;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
57 RZ(a&&w);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
58 d=1; n=1+*AS(a);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
59 v=DAV(a); DO(n-1, d*=*v; v+=n;); d=jfloor(0.5+ABS(d));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
60 if(!d||d>1e20)R w;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
61 v=DAV(w); DO(AN(w), v[i]=jfloor(0.5+d*v[i])/d;);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
62 R w;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
63 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
64
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
65 F1(jtminv){PROLOG;A q,r,*v,y,z;I m,n,*s,t,wr;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
66 F1RANK(2,jtminv,0);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
67 t=AT(w); wr=AR(w); s=AS(w); m=wr?s[0]:1; n=1<wr?s[1]:1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
68 if(!wr)R recip(w);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
69 if(!AN(w)){ASSERT(1==wr||m>=n,EVLENGTH); R cant1(w);}
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
70 if(AN(w)&&t&RAT+XNUM){
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
71 ASSERT(m>=n,EVLENGTH);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
72 if(t&XNUM)RZ(w=cvt(RAT,w));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
73 if(1<wr&&m==n)y=w; else{q=cant1(w); y=pdt(q,w);}
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
74 z=drop(v2(0L,n),gausselm(stitch(y,reshape(v2(n,n),take(sc(1+n),xco1(scf(1.0)))))));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
75 if(2>wr)z=tymes(reshape(mtv,z),w); else if(m>n)z=pdt(z,q);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
76 }else{
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
77 RZ(y=qr(w)); v=AAV(y); q=*v++; r=*v;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
78 z=pdt(rinv(r),t&CMPX?conjug(cant1(q)):cant1(q));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
79 if(t&B01+INT&&2==wr&&m==n)z=icor(r,z);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
80 z=2==wr?z:reshape(shape(w),z);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
81 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
82 EPILOG(z);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
83 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
84
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
85 static B jttridiag(J jt,I n,A a,A x){D*av,d,p,*xv;I i,j,n1=n-1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
86 av=DAV(a); xv=DAV(x); d=xv[0];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
87 for(i=j=0;i<n1;++i){
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
88 ASSERT(d,EVDOMAIN);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
89 p=xv[j+2]/d;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
90 d=xv[j+3]-=p*xv[j+1];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
91 av[i+1]-=p*av[i];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
92 j+=3;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
93 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
94 ASSERT(d,EVDOMAIN);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
95 i=n-1; j=AN(x)-1; av[i]/=d;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
96 for(i=n-2;i>=0;--i){j-=3; av[i]=(av[i]-xv[j+1]*av[i+1])/xv[j];}
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
97 R 1;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
98 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
99
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
100 static F2(jtmdivsp){A a1,x,y;I at,d,m,n,t,*v,xt;P*wp;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
101 ASSERT(2==AR(w),EVRANK);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
102 v=AS(w); n=v[0];
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
103 ASSERT(n>=v[1]&&n==AN(a),EVLENGTH);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
104 ASSERT(n==v[1],EVNONCE);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
105 wp=PAV(w); x=SPA(wp,x); y=SPA(wp,i); a1=SPA(wp,a);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
106 ASSERT(2==AN(a1),EVNONCE);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
107 v=AV(y); m=*AS(y);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
108 ASSERT(m==3*n-2,EVNONCE);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
109 DO(m, d=*v++; d-=*v++; ASSERT(-1<=d&&d<=1,EVNONCE););
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
110 at=AT(a); xt=AT(x); RE(t=maxtype(at,xt)); RE(t=maxtype(t,FL));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
111 RZ(a=cvt(t,a)); RZ(x=cvt(t,x));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
112 if(t&CMPX)RZ(ztridiag(n,a,x)) else RZ(tridiag(n,a,x));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
113 R a;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
114 } /* currently only handles tridiagonal sparse w */
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
115
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
116
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
117 F2(jtmdiv){PROLOG;A q,r,*v,y,z;B b=0;I t;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
118 F2RANK(RMAX,2,jtmdiv,0);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
119 if(AT(a)&SPARSE)RZ(a=denseit(a));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
120 t=AT(w);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
121 if(t&SPARSE)R mdivsp(a,w);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
122 if(t&XNUM+RAT)z=minv(w);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
123 else{
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
124 RZ(y=qr(w)); v=AAV(y); q=*v++; r=*v;
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
125 z=pdt(rinv(r),t&CMPX?conjug(cant1(q)):cant1(q));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
126 b=t&B01+INT&&2==AR(w)&&*AS(w)==*(1+AS(w));
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
127 if(b)z=icor(r,z);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
128 }
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
129 z=pdt(2>AR(w)?reshape(shape(w),z):z,a);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
130 if(b&&AT(a)&B01+INT)z=icor(r,z);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
131 EPILOG(z);
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
diff changeset
132 }