1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include <stdio.h>
8 #include <math.h>
9 #include <stdlib.h>
10 #include <string.h>
11 #include "cs.h"
12 #include "eispack.h"
13 #include "Aomp.h"
14
15 #ifdef isfinite
16 # define IS_GOOD_FLOAT(x) isfinite(x)
17 #else
18 # define IS_GOOD_FLOAT(x) finite(x)
19 # define isfinite finite
20 #endif
21
22 #undef USE_SVDLIB
23 #ifdef USE_SVDLIB
24 #include "svdlib.c"
25 #endif
26
27 void svd_double_ata( int m, int n, double *a, double *s, double *u, double *v ) ;
28
29 /*---------------------------------------------------------------------------*/
30
31 #undef SQR
32 #define SQR(a) ((a)*(a))
33
34 #undef DET3
35 #define DET3(m) ( m[0]*m[4]*m[8]-m[0]*m[7]*m[5]-m[1]*m[3]*m[8] \
36 +m[1]*m[6]*m[5]+m[2]*m[3]*m[7]-m[2]*m[6]*m[4] )
37
38 #undef PI
39 #define PI 3.14159265358979323846
40
41 #undef MIN
42 #define MIN(a,b) (((a)>(b)) ? (b) : (a))
43
44 #undef SWAP
45 #define SWAP(x,y) (th=(x),(x)=(y),(y)=th)
46
47 #undef CSWAP
48 #define CSWAP(i,j) (SWAP(a[i],a[j]),SWAP(a[i+1],a[j+1]),SWAP(a[i+2],a[j+2]))
49
50 #undef EPS
51 #undef EPSQ
52 #define EPS 1.e-8
53 #define EPSQ 1.e-4 /* sqrt(EPS) */
54
55 /*---------------------------------------------------------------------------*/
56 /*! Do a 3x3 symmetric eigen-problem.
57 - INPUT: double a[9] = input matrix; a[i+3*j] = A(i,j) element
58 _ (only the [0], [1], [2], [4], [5], and [8] elements are used)
59 - OUTPUT: e[i] = i'th eigenvalue, with e[0] <= e[1] <= e[2].
60 - OUTPUT: if(dovec) then a[] is replaced with eigenvectors,
61 and this orthogonal matrix will have determinant=1.
62 - METHOD: direct solution of cubic equation
63 - AUTHOR: RW Cox
64 -----------------------------------------------------------------------------*/
65
symeig_3(double * a,double * e,int dovec)66 void symeig_3( double *a , double *e , int dovec )
67 {
68 double aa,bb,cc,dd,ee,ff ;
69 double a1,a2,a3 , qq,rr, qs,th , lam1,lam2,lam3 ;
70 double aba,abb,abc,abd,abe,abf , ann ;
71 double d12,d13,d23 ;
72 double u1,u2,u3 , v1,v2,v3 , w1,w2,w3 , t1,t2,t3 , tn ;
73 double anni ;
74
75 if( a == NULL || e == NULL ) return ;
76
77 /*----- unload matrix into local variables -----*/
78
79 aa = a[0] ; bb = a[1] ; cc = a[2] ; /* matrix is [ aa bb cc ] */
80 dd = a[4] ; ee = a[5] ; ff = a[8] ; /* [ bb dd ee ] */
81 /* [ cc ee ff ] */
82 aba = fabs(aa) ; abb = fabs(bb) ; abc = fabs(cc) ;
83 abd = fabs(dd) ; abe = fabs(ee) ; abf = fabs(ff) ;
84 ann = aba+abb+abc+abd+abe+abf ; /* matrix 'norm' */
85
86 if( ann == 0.0 ){ /* matrix is all zero! */
87 e[0] = e[1] = e[2] = 0.0 ;
88 if( dovec ){
89 a[0] = a[4] = a[8] = 1.0 ;
90 a[1] = a[2] = a[3] = a[5] = a[6] = a[7] = 0.0 ;
91 }
92 return ;
93 }
94
95 /*----- check for matrix that is essentially diagonal -----*/
96
97 if( abb+abc+abe == 0.0 ||
98 ( EPS*aba > (abb+abc) && EPS*abd > (abb+abe) && EPS*abf > (abc+abe) ) ){
99
100 lam1 = aa ; lam2 = dd ; lam3 = ff ;
101
102 if( dovec ){
103 a[0] = a[4] = a[8] = 1.0 ;
104 a[1] = a[2] = a[3] = a[5] = a[6] = a[7] = 0.0 ;
105
106 if( lam1 > lam2 ){ SWAP(lam1,lam2) ; CSWAP(0,3) ; }
107 if( lam1 > lam3 ){ SWAP(lam1,lam3) ; CSWAP(0,6) ; }
108 if( lam2 > lam3 ){ SWAP(lam2,lam3) ; CSWAP(3,6) ; }
109 if( DET3(a) < 0.0 ){ a[6] = -a[6]; a[7] = -a[7]; a[8] = -a[8]; }
110 } else {
111 if( lam1 > lam2 ) SWAP(lam1,lam2) ;
112 if( lam1 > lam3 ) SWAP(lam1,lam3) ;
113 if( lam2 > lam3 ) SWAP(lam2,lam3) ;
114 }
115 e[0] = lam1 ; e[1] = lam2 ; e[2] = lam3 ;
116 return ;
117 }
118
119 /*-- Scale matrix so abs sum is 1; unscale e[i] on output [26 Oct 2005] --*/
120
121 anni = 1.0 / ann ; /* ann != 0, from above */
122 aa *= anni ; bb *= anni ; cc *= anni ;
123 dd *= anni ; ee *= anni ; ff *= anni ;
124
125 /*----- not diagonal ==> must solve cubic polynomial for eigenvalues -----*/
126 /* the cubic polynomial is x**3 + a1*x**2 + a2*x + a3 = 0 */
127
128 a1 = -(aa+dd+ff) ;
129 a2 = (aa*ff+aa*dd+dd*ff - bb*bb-cc*cc-ee*ee) ;
130 a3 = ( aa*(ee*ee-dd*ff) + bb*(bb*ff-cc*ee) + cc*(cc*dd-bb*ee) ) ;
131
132 /*-- Rewrite classical formula for qq as a sum of squares [26 Oct 2005] --*/
133 #if 0
134 qq = (a1*a1 - 3.0*a2) / 9.0 ;
135 #else /* note SQR(a) = (a)*(a) */
136 qq = ( 0.5 * ( SQR(dd-aa) + SQR(ff-aa) + SQR(ff-dd) )
137 + 3.0 * ( bb*bb + cc*cc + ee*ee ) ) / 9.0 ;
138 #endif
139 rr = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3) / 54.0 ;
140
141 if( qq <= 0.0 ){ /*** This should never happen!!! ***/
142 static int nerr=0 ;
143 {
144 if( ++nerr < 4 )
145 fprintf(stderr,"** ERROR in symeig_3: discrim=%g numer=%g\n",qq,rr) ;
146 }
147 qs = qq = rr = 0.0 ;
148 } else {
149 qs = sqrt(qq) ; rr = rr / (qs*qq) ;
150 if( rr < -1.0 ) rr = -1.0 ; else if( rr > 1.0 ) rr = 1.0 ;
151 }
152 th = acos(rr) ;
153
154 lam1 = -2.0 * qs * cos( th /3.0 ) - a1 / 3.0 ;
155 lam2 = -2.0 * qs * cos( (th+2.0*PI)/3.0 ) - a1 / 3.0 ;
156 lam3 = -2.0 * qs * cos( (th+4.0*PI)/3.0 ) - a1 / 3.0 ;
157
158 /*-- if not doing eigenvectors, just sort the eigenvalues to be done --*/
159
160 if( !dovec ){
161 if( lam1 > lam2 ) SWAP(lam1,lam2) ;
162 if( lam1 > lam3 ) SWAP(lam1,lam3) ;
163 if( lam2 > lam3 ) SWAP(lam2,lam3) ;
164 e[0] = ann*lam1 ; e[1] = ann*lam2 ; e[2] = ann*lam3 ; /* re-scale */
165 return ;
166 }
167
168 /*-- are doing eigenvectors; must do double root as a special case --*/
169
170 #undef CROSS /* cross product (x1,x2,x3) X (y1,y2,y3) -> (z1,z2,z3) */
171 #define CROSS(x1,x2,x3,y1,y2,y3,z1,z2,z3) \
172 ( (z1)=(x2)*(y3)-(x3)*(y2), (z2)=(x3)*(y1)-(x1)*(y3), (z3)=(x1)*(y2)-(x2)*(y1) )
173
174 d12 = fabs(lam1-lam2) ; d13 = fabs(lam1-lam3) ; d23 = fabs(lam2-lam3) ;
175 rr = MIN(d12,d13) ; rr = MIN(rr,d23) ;
176
177 if( rr > EPS*ann ){ /*---- not a double root ----*/
178
179 if( lam1 > lam2 ) SWAP(lam1,lam2) ; /* start by sorting eigenvalues */
180 if( lam1 > lam3 ) SWAP(lam1,lam3) ;
181 if( lam2 > lam3 ) SWAP(lam2,lam3) ;
182 e[0] = ann*lam1 ; e[1] = ann*lam2 ; e[2] = ann*lam3 ; /* re-scale */
183
184 /* find eigenvector for lam1 by computing Ay-lam1*y for
185 vectors y1=[1,0,0], y2=[0,1,0], and y3=[0,0,1]; the eigenvector
186 is orthogonal to all of these, so use the cross product to get it */
187
188 u1 = aa-lam1 ; u2 = bb ; u3 = cc ; /* A*y1 - lam1*y1 */
189 v1 = bb ; v2 = dd-lam1 ; v3 = ee ; /* A*y2 - lam1*y2 */
190 CROSS(u1,u2,u3 , v1,v2,v3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
191 if( tn < EPSQ*ann ){ /* u and v were parallel? */
192 w1 = cc ; w2 = ee ; w3 = ff-lam1 ; /* A*y3 - lam1*y3 */
193 CROSS(u1,u2,u3 , w1,w2,w3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
194 if( tn < EPSQ*ann ){ /* u and w were parallel? */
195 CROSS(v1,v2,v3 , w1,w2,w3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
196 }
197 }
198 a[0] = t1/tn ; a[1] = t2/tn ; a[2] = t3/tn ; /* normalize */
199
200 /* do same for lam2 */
201
202 u1 = aa-lam2 ; u2 = bb ; u3 = cc ;
203 v1 = bb ; v2 = dd-lam2 ; v3 = ee ;
204 CROSS(u1,u2,u3 , v1,v2,v3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
205 if( tn < EPSQ*ann ){
206 w1 = cc ; w2 = ee ; w3 = ff-lam2 ;
207 CROSS(u1,u2,u3 , w1,w2,w3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
208 if( tn < EPSQ*ann ){
209 CROSS(v1,v2,v3 , w1,w2,w3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
210 }
211 }
212 a[3] = t1/tn ; a[4] = t2/tn ; a[5] = t3/tn ; /* normalize */
213
214 /* orthgonality of eigenvectors ==> can get last one by cross product */
215
216 #if 1
217 CROSS( a[0],a[1],a[2] , a[3],a[4],a[5] , a[6],a[7],a[8] ) ;
218 #else
219 u1 = aa-lam3 ; u2 = bb ; u3 = cc ;
220 v1 = bb ; v2 = dd-lam3 ; v3 = ee ;
221 CROSS(u1,u2,u3 , v1,v2,v3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
222 if( tn < EPSQ*ann ){
223 w1 = cc ; w2 = ee ; w3 = ff-lam3 ;
224 CROSS(u1,u2,u3 , w1,w2,w3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
225 if( tn < EPSQ*ann ){
226 CROSS(v1,v2,v3 , w1,w2,w3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
227 }
228 }
229 a[6] = t1/tn ; a[7] = t2/tn ; a[8] = t3/tn ; /* normalize */
230 #endif
231
232 } else { /*---- if here, we have a double root ----*/
233
234 /* make sure that we have lam1=lam2 and lam3 is the outlier */
235
236 if( d13 < d12 && d13 < d23 ) SWAP(lam2,lam3) ;
237 else if( d23 < d12 && d23 < d13 ) SWAP(lam1,lam3) ;
238 lam1 = lam2 = 0.5*(lam1+lam2) ;
239
240 /* compute eigenvector for lam3 using method as above */
241
242 u1 = aa-lam3 ; u2 = bb ; u3 = cc ;
243 v1 = bb ; v2 = dd-lam3 ; v3 = ee ;
244 CROSS(u1,u2,u3 , v1,v2,v3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
245 if( tn < EPSQ*ann ){
246 w1 = cc ; w2 = ee ; w3 = ff-lam3 ;
247 CROSS(u1,u2,u3 , w1,w2,w3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
248 if( tn < EPSQ*ann ){
249 CROSS(v1,v2,v3 , w1,w2,w3 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
250 }
251 }
252 w1 = a[6] = t1/tn ; w2 = a[7] = t2/tn ; w3 = a[8] = t3/tn ;
253
254 /* find a vector orthogonal to it */
255
256 CROSS(w1,w2,w3 , 1,0,0 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
257 if( tn < EPSQ ){
258 CROSS(w1,w2,w3 , 0,1,0 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
259 if( tn < EPSQ ){
260 CROSS(w1,w2,w3 , 0,0,1 , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
261 }
262 }
263 a[0] = t1/tn ; a[1] = t2/tn ; a[2] = t3/tn ;
264
265 /* and the final vector is the cross product of these two */
266
267 CROSS( w1,w2,w3 , a[0],a[1],a[2] , t1,t2,t3 ) ; tn = sqrt(t1*t1+t2*t2+t3*t3) ;
268 a[3] = t1/tn ; a[4] = t2/tn ; a[5] = t3/tn ;
269
270 /* sort results (we know lam1==lam2) */
271
272 if( lam1 > lam3 ){
273 SWAP(lam1,lam3) ; CSWAP(0,6) ;
274 if( DET3(a) < 0.0 ){ a[6] = -a[6]; a[7] = -a[7]; a[8] = -a[8]; }
275 }
276
277 e[0] = ann*lam1 ; e[1] = ann*lam2 ; e[2] = ann*lam3 ; /* re-scale */
278 }
279
280 return ;
281 }
282
283 /*---------------------------------------------------------------------------*/
284 /*! 2x2 symmetric eigenvalue/vector problem, like symeig_3() above. */
285
symeig_2(double * a,double * e,int dovec)286 void symeig_2( double *a , double *e , int dovec )
287 {
288 double sxx,sxy,syy , lam1,lam2 , ss,tt , x,y ;
289
290 if( a == NULL || e == NULL ) return ;
291
292 /*----- unload matrix into local variables -----*/
293
294 sxx = a[0] ; sxy = a[1] ; syy = a[3] ;
295
296 ss = fabs(sxx) ; tt = fabs(syy) ; if( ss > tt ) ss = tt ;
297
298 if( fabs(sxy) < EPS*ss ){ /*--- essentially a diagonal matrix ---*/
299 if( sxx <= syy ){
300 lam1 = sxx ; lam2 = syy ;
301 if( dovec ){ a[0]=a[3]=1.0; a[1]=a[2]=0.0; }
302 } else {
303 lam1 = syy ; lam2 = sxx ;
304 if( dovec ){ a[0]=a[3]=1.0 ; a[1]=a[2]=1.0; }
305 }
306 e[0] = lam1 ; e[1] = lam2 ;
307 return ;
308 }
309
310 /*--- non-diagonal matrix ==> solve quadratic equation for eigenvalues ---*/
311
312 ss = sqrt( (sxx-syy)*(sxx-syy) + 4.0*sxy*sxy ) ; /* positive */
313 lam1 = 0.5 * ( sxx + syy - ss ) ; /* smaller */
314 lam2 = 0.5 * ( sxx + syy + ss ) ; /* larger */
315
316 if( dovec ){
317 x = 2.0*sxy ; y = syy - sxx - ss ; tt = sqrt(x*x+y*y) ;
318 a[0] = x/tt ; a[1] = y/tt ;
319
320 y = syy - sxx + ss ; tt = sqrt(x*x+y*y) ;
321 a[2] = x/tt ; a[3] = y/tt ;
322 }
323 e[0] = lam1 ; e[1] = lam2 ;
324 return ;
325 }
326
327 /*--------------------------------------------------------------------------*/
328
329 static int forbid_23 = 0 ;
330
331 /*! To turn off use of symeig_3() and symeig_2() in the symeig functions. */
332
symeig_forbid_23(int ii)333 void symeig_forbid_23( int ii ){ forbid_23 = ii ; return ; }
334
335 /*--------------------------------------------------------------------------
336 Compute the eigenvalue/vector decomposition of a symmetric matrix,
337 stored in double precision.
338 n = order of matrix
339 a = on input: matrix(i,j) is in a[i+n*j] for i=0..n-1 , j=0..n-1
340 output: a[i+n*j] has the i'th component of the j'th eigenvector
341 e = on input: not used (but the calling program must
342 allocate the space for e[0..n-1])
343 output: e[j] has the j'th eigenvalue, ordered so that
344 e[0] <= e[1] <= ... <= e[n-1]
345 Uses the f2c translation of EISPACK.
346 ----------------------------------------------------------------------------*/
347
symeig_double(int n,double * a,double * e)348 void symeig_double( int n , double *a , double *e )
349 {
350 integer nm , matz , ierr ;
351 double *fv1 , *fv2 ;
352
353 if( a == NULL || e == NULL || n < 1 ) return ;
354
355 /* special cases of small n (much faster than EISPACK) */
356
357 if( n == 1 ){
358 e[0] = a[0] ; a[0] = 1.0 ; return ; /* degenerate case */
359 } else if( !forbid_23 ){
360 if( n == 2 ){ symeig_2( a , e , 1 ) ; return ; }
361 if( n == 3 ){ symeig_3( a , e , 1 ) ; return ; }
362 }
363
364 /*-- default code: n > 3 --*/
365
366 fv1 = (double *) malloc(sizeof(double)*n) ; /* workspaces */
367 fv2 = (double *) malloc(sizeof(double)*n) ;
368
369 nm = n ; matz = 1 ; ierr = 0 ;
370
371 rs_( &nm , &nm , a , e , &matz , a , fv1 , fv2 , &ierr ) ;
372
373 free((void *)fv1) ; free((void *)fv2) ;
374 return ;
375 }
376
377 /*------------------ just compute the eigenvalues -------------------*/
378
symeigval_double(int n,double * a,double * e)379 void symeigval_double( int n , double *a , double *e )
380 {
381 integer nm , matz , ierr ;
382 double *fv1 , *fv2 ;
383
384 if( a == NULL || e == NULL || n < 1 ) return ;
385
386 /* special cases of small n (much faster than EISPACK) */
387
388 if( n == 1 ){
389 e[0] = a[0] ; return ; /* degenerate case */
390 } else if( !forbid_23 ){
391 if( n == 2 ){ symeig_2( a , e , 0 ) ; return ; }
392 if( n == 3 ){ symeig_3( a , e , 0 ) ; return ; }
393 }
394
395 /*-- here, deal with general n > 3 --*/
396
397 fv1 = (double *) malloc(sizeof(double)*(n+9)) ; /* workspaces */
398 fv2 = (double *) malloc(sizeof(double)*(n+9)) ;
399
400 nm = n ; matz = 0 ; ierr = 0 ;
401
402 rs_( &nm , &nm , a , e , &matz , a , fv1 , fv2 , &ierr ) ;
403
404 if( ierr != 0 )
405 fprintf(stderr,"** ERROR: symeigval_double error code = %d\n",(int)ierr) ;
406
407 free((void *)fv1) ; free((void *)fv2) ;
408 return ;
409 }
410
411 /*--------------------------------------------------------------------*/
412 /*! Return eigenvalues/eigenvectors indexed from bb to tt (bb <= tt),
413 where index #0 = smallest eigenvalue, index #n-1 = largest
414 n = order of matrix
415 a = on input: matrix(i,j) is in a[i+n*j] for i=0..n-1 , j=0..n-1
416 output: a[i+n*j] has the i'th component of the j'th
417 eigenvector, for j=0..tt-bb.
418 However, if novec!=0, then eigenvectors will not be computed;
419 only the eigenvalues will be output.
420 e = on input: not used (but the calling program must
421 allocate the space for e[0..tt-bb])
422 output: e[j] has the j'th eigenvalue, ordered so that
423 e[0] <= e[1] <= ... <= e[tt-bb]
424
425 Return value is 0 for all being done OK, nonzero for error.
426 *//*------------------------------------------------------------------*/
427
symeig_irange(int n,double * a,double * e,int bb,int tt,int novec)428 int symeig_irange( int n, double *a, double *e, int bb, int tt, int novec )
429 {
430 integer nm , m11,mmm , ierr , *ind ;
431 double *fv1, *fv2, *fv3, eps1, lb,ub, *rv4,*rv5,*rv6,*rv7,*rv8, *zzz ;
432 int ii , nval ;
433
434 if( n < 1 || a == NULL || e == NULL || bb < 0 || tt < bb || tt >= n )
435 return -66666 ;
436
437 if( bb==0 && tt==n-1 ){ symeig_double( n , a , e ) ; return 0 ; }
438
439 /* reduction to tridiagonal form (stored in fv1..3) */
440
441 nm = n ;
442 fv1 = (double *) malloc(sizeof(double)*(n+9)) ; /* workspaces */
443 fv2 = (double *) malloc(sizeof(double)*(n+9)) ;
444 fv3 = (double *) malloc(sizeof(double)*(n+9)) ;
445
446 tred1_( &nm , &nm , a , fv1,fv2,fv3 ) ;
447
448 /* determination of the desired eigenvalues of the tridiagonal matrix */
449
450 eps1 = 0.0 ;
451 m11 = bb+1 ;
452 mmm = tt-bb+1 ;
453 ierr = 0 ;
454 ind = (integer *)malloc(sizeof(integer)*(n+9)) ;
455 rv4 = (double *) malloc(sizeof(double) *(n+9)) ;
456 rv5 = (double *) malloc(sizeof(double) *(n+9)) ;
457
458 tridib_( &nm , &eps1 , fv1,fv2,fv3 , &lb,&ub , &m11,&mmm , e ,
459 ind , &ierr , rv4,rv5 ) ;
460
461 if( ierr != 0 || novec != 0 ){
462 free(rv5); free(rv4); free(ind); free(fv3); free(fv2); free(fv1);
463 return -ierr ;
464 }
465
466 /* determination of the eigenvectors of the tridiagonal matrix */
467
468 nval = nm * mmm ;
469 zzz = (double *) malloc(sizeof(double) *nval ) ;
470 rv6 = (double *) malloc(sizeof(double) *(n+9)) ;
471 rv7 = (double *) malloc(sizeof(double) *(n+9)) ;
472 rv8 = (double *) malloc(sizeof(double) *(n+9)) ;
473
474 tinvit_( &nm , &nm , fv1,fv2,fv3 , &mmm , e ,
475 ind , zzz , &ierr , rv4,rv5,rv6,rv7,rv8 ) ;
476
477 if( ierr != 0 ){
478 free(rv8); free(rv7); free(rv6); free(zzz);
479 free(rv5); free(rv4); free(ind); free(fv3); free(fv2); free(fv1);
480 return ierr ;
481 }
482
483 /* transform eigenvectors back to original space */
484
485 trbak1_( &nm , &nm , a , fv2 , &mmm , zzz ) ;
486
487 /* copy output eigenvectors into a */
488
489 for( ii=0 ; ii < nval ; ii++ ) a[ii] = zzz[ii] ;
490
491 free(rv8); free(rv7); free(rv6); free(zzz);
492 free(rv5); free(rv4); free(ind); free(fv3); free(fv2); free(fv1);
493 return 0 ;
494 }
495
496 /*----------------------------------------------------------------------------*/
497
498 #undef A
499 #define A(i,j) asym[(i)+(j)*nsym]
500
501 /*----------------------------------------------------------------------------*/
502 /*! Compute the nvec principal singular vectors of a set of m columns, each
503 of length n, stored in array xx[i+j*n] for i=0..n-1, j=0..m-1.
504 The singular values (largest to smallest) are stored in sval, and
505 the left singular vectors [first nvec columns of U in X = U S V'] are
506 stored into uvec[i+j*n] for i=0..n-1, j=0..nvec-1.
507
508 The return value is the number of vectors computed. If the return
509 value is not positive, something bad happened. Normally, the return
510 value would be the same as nvec, but it cannot be larger than MIN(n,m).
511
512 If sval==NULL, then the output into sval is skipped.
513 If uval==NULL, then the output into uval is skipped.
514 If both are NULL, exactly why did you want to call this function?
515 *//*--------------------------------------------------------------------------*/
516
first_principal_vectors(int n,int m,float * xx,int nvec,float * sval,float * uvec)517 int first_principal_vectors( int n , int m , float *xx ,
518 int nvec , float *sval , float *uvec )
519 {
520 int nn=n , mm=m , nsym , ii,jj,kk,qq ;
521 double *asym , *deval ;
522 register double sum , qsum ; register float *xj , *xk ;
523
524 nsym = MIN(nn,mm) ; /* size of the symmetric matrix to create */
525
526 if( nsym < 1 || xx == NULL || (uvec == NULL && sval == NULL) ) return -666 ;
527
528 if( nvec > nsym ) nvec = nsym ; /* can't compute more vectors than nsym! */
529
530 { asym = (double *)malloc(sizeof(double)*nsym*nsym) ; /* symmetric matrix */
531 deval = (double *)malloc(sizeof(double)*nsym) ; /* its eigenvalues */
532 }
533
534 /** setup matrix to eigensolve: choose smaller of [X]'[X] and [X][X]' **/
535 /** since [X] is n x m, [X]'[X] is m x m and [X][X]' is n x n **/
536
537 if( nn > mm ){ /* more rows than columns: */
538 /* so [A] = [X]'[X] = m x m */
539 int n1 = nn-1 ;
540 for( jj=0 ; jj < mm ; jj++ ){
541 xj = xx + jj*nn ;
542 for( kk=0 ; kk <= jj ; kk++ ){
543 sum = 0.0 ; xk = xx + kk*nn ;
544 for( ii=0 ; ii < n1 ; ii+=2 ) sum += xj[ii]*xk[ii] + xj[ii+1]*xk[ii+1];
545 if( ii == n1 ) sum += xj[ii]*xk[ii] ;
546 A(jj,kk) = sum ; if( kk < jj ) A(kk,jj) = sum ;
547 }
548 }
549
550 } else { /* more columns than rows: */
551 /* so [A] = [X][X]' = n x n */
552 float *xt ; int m1=mm-1 ;
553 xt = (float *)malloc(sizeof(float)*nn*mm) ;
554 for( jj=0 ; jj < mm ; jj++ ){ /* form [X]' into array xt */
555 for( ii=0 ; ii < nn ; ii++ ) xt[jj+ii*mm] = xx[ii+jj*nn] ;
556 }
557
558 for( jj=0 ; jj < nn ; jj++ ){
559 xj = xt + jj*mm ;
560 for( kk=0 ; kk <= jj ; kk++ ){
561 sum = 0.0 ; xk = xt + kk*mm ;
562 for( ii=0 ; ii < m1 ; ii+=2 ) sum += xj[ii]*xk[ii] + xj[ii+1]*xk[ii+1];
563 if( ii == m1 ) sum += xj[ii]*xk[ii] ;
564 A(jj,kk) = sum ; if( kk < jj ) A(kk,jj) = sum ;
565 }
566 }
567
568 free(xt) ; /* don't need this no more */
569 }
570
571 /** compute the nvec eigenvectors corresponding to largest eigenvalues **/
572 /** these eigenvectors are stored on top of first nvec columns of asym **/
573
574 ii = symeig_irange( nsym, asym, deval, nsym-nvec, nsym-1, (uvec==NULL) ) ;
575
576 if( ii != 0 ){
577 { free(deval) ; free(asym) ; }
578 return -33333 ; /* eigensolver failed!? */
579 }
580
581 /** Store singular values (sqrt of eigenvalues), if desired:
582 Note that symeig_irange returns things smallest to largest,
583 but we want largest to smallest, so have to reverse the order **/
584
585 if( sval != NULL ){
586 for( jj=0 ; jj < nvec ; jj++ ){
587 sum = deval[nvec-1-jj] ;
588 sval[jj] = (sum <= 0.0) ? 0.0 : sqrt(sum) ;
589 }
590 }
591
592 /** if no output vectors desired, we are done done done!!! **/
593
594 if( uvec == NULL ){
595 { free(deval) ; free(asym) ; }
596 return nvec ;
597 }
598
599 /** SVD is [X] = [U] [S] [V]', where [U] = desired output vectors
600
601 case n <= m: [A] = [X][X]' = [U] [S][S]' [U]'
602 so [A][U] = [U] [S][S]'
603 so eigenvectors of [A] are just [U]
604
605 case n > m: [A] = [X]'[X] = [V] [S]'[S] [V]'
606 so [A][V] = [V] [S'][S]
607 so eigenvectors of [A] are [V], but we want [U]
608 note that [X][V] = [U] [S]
609 so pre-multiplying each column vector in [V] by matrix [X]
610 will give the corresponding column in [U], but scaled;
611 below, just L2-normalize the column to get output vector **/
612
613 if( nn <= mm ){ /* copy eigenvectors into output directly */
614 /* (e.g., more vectors than time points) */
615 for( jj=0 ; jj < nvec ; jj++ ){
616 qq = nvec-1-jj ; /* eigenvalues are in reversed order */
617 for( ii=0 ; ii < nn ; ii++ )
618 uvec[ii+jj*nn] = (float)asym[ii+qq*nn] ;
619 }
620
621 } else { /* n > m: transform eigenvectors to get left singular vectors */
622 /* (e.g., more time points than vectors) */
623
624 for( jj=0 ; jj < nvec ; jj++ ){
625 qq = nvec-1-jj ; qsum = 0.0 ; /* eigenvalues are in reversed order */
626 for( ii=0 ; ii < nn ; ii++ ){
627 sum = 0.0 ;
628 for( kk=0 ; kk < mm ; kk++ ) sum += xx[ii+kk*nn] * asym[kk+qq*mm] ;
629 uvec[ii+jj*nn] = sum ; qsum += sum*sum ;
630 }
631 if( qsum > 0.0 ){ /* L2 normalize */
632 register float fac ;
633 fac = (float)(1.0/sqrt(qsum)) ;
634 for( ii=0 ; ii < nn ; ii++ ) uvec[ii+jj*nn] *= fac ;
635 }
636 }
637 }
638
639 /** free at last!!! **/
640
641 { free(deval) ; free(asym) ; }
642 return nvec ;
643 }
644
645 #undef A
646
647 /*--------------------------------------------------------------------*/
648
649 #define CHECK_SVD /* unfortunately, this is necessary */
650
651 #undef CHK
652 #ifdef CHECK_SVD
653 # define CHK 1
654 # define A(i,j) aa[(i)+(j)*m]
655 # define U(i,j) uu[(i)+(j)*m]
656 # define V(i,j) vv[(i)+(j)*n]
657 #else
658 # define CHK 0
659 #endif
660
661 /** setup for sorting SVD values:
662 0 = no sort (whatever the function returns)
663 +1 = sort in increasing order of singular values
664 -1 = sort in descending order of singular values **/
665
666 static int svd_sort = 0 ;
set_svd_sort(int ss)667 void set_svd_sort( int ss ){ svd_sort = ss; }
668
669 /*----------------------------------------------------------------------------*/
670 /*! Compute SVD of double precision matrix: T
671 [a] = [u] diag[s] [v]
672 - m = # of rows in a = length of each column
673 - n = # of columns in a = length of each row
674 - a = pointer to input matrix; a[i+j*m] has the (i,j) element
675 (m X n matrix, stored in column-first order)
676 - s = pointer to output singular values; length = n (cannot be NULL)
677 - u = pointer to output matrix, if desired; length = m*n (m X n matrix)
678 - v = pointer to output matrix, if desired; length = n*n (n x n matrix)
679
680 Modified 10 Jan 2007 to add sorting of s and corresponding columns of u & v.
681 ------------------------------------------------------------------------------*/
682
svd_double(int m,int n,double * a,double * s,double * u,double * v)683 void svd_double( int m, int n, double *a, double *s, double *u, double *v )
684 {
685 integer mm,nn , lda,ldu,ldv , ierr ;
686 doublereal *aa, *ww , *uu , *vv , *rv1 ;
687 logical matu , matv ;
688
689 if( a == NULL || s == NULL || m < 1 || n < 1 ) return ;
690
691 mm = m ;
692 nn = n ;
693 aa = a ;
694 lda = m ;
695 ww = s ;
696
697 /* make space for u matrix, if not supplied */
698
699 if( u == NULL ){
700 matu = (logical) CHK ;
701 uu = (doublereal *)calloc(sizeof(double),m*n) ;
702 } else {
703 matu = (logical) 1 ;
704 uu = u ;
705 }
706 ldu = m ;
707
708 /* make space for v matrix if not supplied */
709
710 if( v == NULL ){
711 matv = (logical) CHK ;
712 vv = (CHK) ? (doublereal *)calloc(sizeof(double),n*n) : NULL ;
713 } else {
714 matv = (logical) 1 ;
715 vv = v ;
716 }
717 ldv = n ;
718
719 rv1 = (double *)calloc(sizeof(double),n) ; /* workspace */
720
721 /** the actual SVD **/
722
723 (void) svd_( &mm , &nn , &lda , aa , ww ,
724 &matu , &ldu , uu , &matv , &ldv , vv , &ierr , rv1 ) ;
725
726 /** the following code is to check if the SVD worked **/
727
728 #ifdef CHECK_SVD
729 /** back-compute [A] from [U] diag[ww] [V]'
730 and see if it is close to the input matrix;
731 if not, compute the results in another function;
732 this is needed because the svd() function compiles with
733 rare computational errors on some compilers' optimizers **/
734 { register int i,j,k ; register doublereal aij ; double err=0.0,amag=1.e-12 ;
735 for( j=0 ; j < n ; j++ ){
736 for( i=0 ; i < m ; i++ ){
737 aij = A(i,j) ; amag += fabs(aij) ;
738 for( k=0 ; k < n ; k++ ) aij -= U(i,k)*V(j,k)*ww[k] ;
739 err += fabs(aij) ;
740 }}
741 amag /= (m*n) ; /* average absolute value of matrix elements */
742 err /= (m*n) ; /* average absolute error per matrix element */
743 if( err >= 1.e-5*amag || !IS_GOOD_FLOAT(err) ){
744 fprintf(stderr,"\n **** SVD avg err=%g; recomputing ...",err) ;
745
746 #if 1 /* perturb matrix slightly */
747 { double arep=1.e-13*amag , *aj ;
748 for( j=0 ; j < nn ; j++ ){
749 aj = aa + j*mm ;
750 #if 1
751 for( i=0 ; i < mm ; i++ ) aj[i] += (drand48()-0.5)*arep ;
752 #else
753 for( i=0 ; i < mm ; i++ ) if( aj[i] != 0.0 ) break ;
754 if( i == mm ){
755 for( i=0 ; i < mm ; i++ ) aj[i] = (drand48()-0.5)*arep ;
756 }
757 #endif
758 }
759 }
760 #endif
761
762 /* svd_slow is compiled without optimization */
763
764 (void) svd_slow_( &mm , &nn , &lda , aa , ww ,
765 &matu , &ldu , uu , &matv , &ldv , vv , &ierr , rv1 ) ;
766 err = 0.0 ;
767 for( j=0 ; j < n ; j++ ){
768 for( i=0 ; i < m ; i++ ){
769 aij = A(i,j) ;
770 for( k=0 ; k < n ; k++ ) aij -= U(i,k)*V(j,k)*ww[k] ;
771 err += fabs(aij) ;
772 }}
773 err /= (m*n) ;
774
775 /* not fixed YET? try another algorithm (one that's usually slower) */
776
777 if( err >= 1.e-5*amag || !IS_GOOD_FLOAT(err) ){
778 fprintf(stderr," new avg err=%g; re-recomputing the hard way ...",err) ;
779 #ifdef USE_SVDLIB
780 fprintf(stderr,"\n") ;
781 SVDVerbosity = 2 ;
782 AFNI_svdLAS2( mm , nn , aa , ww , uu , vv ) ; /* svdlib.c */
783 SVDVerbosity = 0 ;
784 #else
785 svd_double_ata( mm , nn , aa , ww , uu, vv ) ; /* below */
786 #endif
787 err = 0.0 ;
788 for( j=0 ; j < n ; j++ ){
789 for( i=0 ; i < m ; i++ ){
790 aij = A(i,j) ;
791 for( k=0 ; k < n ; k++ ) aij -= U(i,k)*V(j,k)*ww[k] ;
792 err += fabs(aij) ;
793 }}
794 err /= (m*n) ;
795 fprintf(stderr," newer avg err=%g %s" ,
796 err ,
797 (err >= 1.e-5*amag || !IS_GOOD_FLOAT(err)) ? "**BAD** :-(" : "**OK** :-)" ) ;
798 } else {
799 fprintf(stderr," new avg error=%g **OK** :-)",err) ;
800 }
801 fprintf(stderr,"\n\n") ;
802 }
803 }
804 #endif
805
806 free((void *)rv1) ;
807
808 /* discard [u] and [v] spaces if not needed for output */
809
810 if( u == NULL && uu != NULL ) free((void *)uu) ;
811 if( v == NULL && vv != NULL ) free((void *)vv) ;
812
813 /*--- 10 Jan 2007: sort the singular values and columns of U and V ---*/
814
815 if( n > 1 && svd_sort != 0 ){
816 double *sv , *uv ; int *iv , jj,kk ;
817 sv = (double *)malloc(sizeof(double)*n) ;
818 iv = (int *) malloc(sizeof(int) *n) ;
819 for( kk=0 ; kk < n ; kk++ ){
820 iv[kk] = kk ; sv[kk] = (svd_sort > 0) ? s[kk] : -s[kk] ;
821 }
822 qsort_doubleint( n , sv , iv ) ;
823 if( u != NULL ){
824 double *cc = (double *)calloc(sizeof(double),m*n) ;
825 #pragma omp critical
826 (void)AAmemcpy( cc , u , sizeof(double)*m*n ) ;
827 for( jj=0 ; jj < n ; jj++ ){
828 kk = iv[jj] ; /* where the new jj-th col came from */
829 #pragma omp critical
830 (void)AAmemcpy( u+jj*m , cc+kk*m , sizeof(double)*m ) ;
831 }
832 free((void *)cc) ;
833 }
834 if( v != NULL ){
835 double *cc = (double *)calloc(sizeof(double),n*n) ;
836 #pragma omp critical
837 (void)AAmemcpy( cc , v , sizeof(double)*n*n ) ;
838 for( jj=0 ; jj < n ; jj++ ){
839 kk = iv[jj] ;
840 #pragma omp critical
841 (void)AAmemcpy( v+jj*n , cc+kk*n , sizeof(double)*n ) ;
842 }
843 free((void *)cc) ;
844 }
845 for( kk=0 ; kk < n ; kk++ )
846 s[kk] = (svd_sort > 0) ? sv[kk] : -sv[kk] ;
847 free((void *)iv) ; free((void *)sv) ;
848 }
849
850 return ;
851 }
852
853 /*----------------------------------------------------------------------------*/
854 /* SVD via eigensolution of A'A (presumes n < m) */
855
svd_double_ata(int m,int n,double * a,double * s,double * u,double * v)856 void svd_double_ata( int m, int n, double *a, double *s, double *u, double *v )
857 {
858 double *ata , *ai,*aj , sum ;
859 int ii,jj,kk ;
860
861 if( a == NULL || s == NULL || m < 1 || n < 1 ) return ;
862
863 /* make A'A matrix */
864
865 ata = (double *)malloc(sizeof(double)*n*n) ;
866 for( ii=0 ; ii < n ; ii++ ){
867 ai = a + ii*m ; /* ii-th column of A */
868 for( jj=0 ; jj <= ii ; jj++ ){
869 aj = a + jj*m ; /* jj-th column of A */
870 for( sum=0.0,kk=0 ; kk < m ; kk++ ) sum += ai[kk]*aj[kk] ;
871 ata[ii+jj*n] = sum ; if( jj < ii ) ata[jj+ii*n] = sum ;
872 }
873 }
874
875 /* Since A = U S V', then A'A = V S^2 V', or A'A V = V S^2,
876 so the eigensolution of A'A gives S^2 as the eigenvalues,
877 and gives V as the eigenvectors (which here over-write ata). */
878
879 symeig_double( n , ata , s ) ;
880
881 /* sqrt the eigenvalues to get the singular values */
882
883 for( ii=0 ; ii < n ; ii++ )
884 s[ii] = (s[ii] <= 0.0) ? 0.0 : sqrt(s[ii]) ;
885
886 /* copy ata into the output place for V */
887
888 if( v != NULL ){
889 #pragma omp critical
890 AAmemcpy( v , ata , sizeof(double)*n*n ) ;
891 }
892
893 /* manufacture U from A and V (which is in ata, recall):
894 AV = US, so if we form the matrix product AV, then
895 L2-normalizing each column will give us the matrix U. */
896
897 if( u != NULL ){
898 double *uj , *vj ;
899 for( jj=0 ; jj < n ; jj++ ){
900 uj = u + jj*m ; /* jj-th column of U */
901 vj = ata + jj*n ; /* jj-th column of V */
902 for( ii=0 ; ii < m ; ii++ ){ /* A_{ii,kk} V_{kk,jj} */
903 for( sum=0.0,kk=0 ; kk < n ; kk++ ) sum += a[ii+kk*m]*vj[kk] ;
904 uj[ii] = sum ;
905 }
906 for( sum=0.0,ii=0 ; ii < m ; ii++ ) sum += uj[ii]*uj[ii] ;
907 if( sum > 0.0 ){
908 for( sum=1.0/sqrt(sum),ii=0 ; ii < m ; ii++ ) uj[ii] *= sum ;
909 }
910 }
911 }
912
913 free(ata) ; return ;
914 }
915
916 /*--------------------------------------------------------------------*/
917
918 #define DBG_PC 0
919 /*!
920 Calculate the covariance matrix of data_mat based on code in 3dpc
921
922 data_mat: Data matrix containg num_cols vectors that have num_rows elements.
923 cov_mat: On output, cov_mat will contain the covariance matrix. Caller must
924 allocate num_cols x num_cols elements for it.
925 Both matrices are stored in column major order. M(i,j) = Mv(i+j*num_rows);
926 row_mask: mask vector num_rows long.
927 NULL for no masking.
928 num_rows, num_cols: 1st and 2nd dimensions of data_mat
929 norm: flag for normalizing covariance.
930 -1 no normalization,
931 0, normalize by num_rows - 1,
932 1, normalize by num_rows
933 remove_mean: 0 : do nothing
934 1 : remove the mean of each column in data_mat (like -dmean in 3dpc)
935
936 To match matlab's cov function, you need to remove the mean of each column
937 and set norm to 0 (default for matlab) or 1
938
939 the function returns the trace of the covariance matrix if all went well;
940 returns -1.0 in case of error.
941
942 */
943
covariance(float * data_mat,double * cov_mat,unsigned char * row_mask,int num_rows,int num_cols,int norm,int remove_mean,int be_quiet)944 double covariance(float *data_mat, double *cov_mat, unsigned char * row_mask,
945 int num_rows, int num_cols, int norm, int remove_mean, int be_quiet)
946 {
947 double atrace, dsum, normval=0.0;
948 int idel, jj, nn, mm, ifirst, ilast, ii, PC_be_quiet, kk, nsum;
949
950 if (norm == 0) normval = (double)num_rows - 1.0;
951 else if (norm == 1) normval = (double)num_rows;
952 else if (norm == -1) normval = 0.0f;
953 else {
954 fprintf(stderr,"*** norm value of %d is not acceptable.\n", norm); return(-1.0);
955 }
956
957 if (remove_mean == 1) {
958 for( jj=0 ; jj < num_cols ; jj++ ){ /* for each column */
959 nsum = 0;
960 dsum = 0.0 ;
961 if( row_mask == NULL ){
962 for( kk=0 ; kk < num_rows ; kk++ ) dsum += data_mat[kk+jj*num_rows] ;
963 dsum /= (double)num_rows;
964 } else {
965 for( kk=0 ; kk < num_rows ; kk++ )
966 if( row_mask[kk] ) { dsum += data_mat[kk+jj*num_rows] ; ++nsum; }
967 dsum /= (double)nsum;
968 }
969 if( row_mask == NULL ){
970 for( kk=0 ; kk < num_rows ; kk++ ) data_mat[kk+jj*num_rows] -= dsum ;
971 } else {
972 for( kk=0 ; kk < num_rows ; kk++ )
973 if( row_mask[kk] ) { data_mat[kk+jj*num_rows] -= dsum; }
974 }
975 }
976 }
977
978 idel = 1 ; /* ii goes forward */
979 for( jj=0 ; jj < num_cols ; jj++ ){
980
981 ifirst = (idel==1) ? 0 : jj ; /* back and forth in ii to */
982 ilast = (idel==1) ? jj+1 : -1 ; /* maximize use of cache/RAM */
983
984 for( ii=ifirst ; ii != ilast ; ii += idel ){
985 dsum = 0.0 ;
986 if( row_mask == NULL ){
987 for( kk=0 ; kk < num_rows ; kk++ ) dsum += data_mat[kk+ii*num_rows] * data_mat[kk+jj*num_rows] ;
988 } else {
989 for( kk=0 ; kk < num_rows ; kk++ )
990 if( row_mask[kk] ) dsum += data_mat[kk+ii*num_rows] * data_mat[kk+jj*num_rows] ;
991 }
992 if (normval > 1) cov_mat[ii+jj*num_cols] = cov_mat[jj+ii*num_cols] = dsum/normval ;
993 else cov_mat[ii+jj*num_cols] = cov_mat[jj+ii*num_cols] = dsum;
994 }
995
996 if( !be_quiet ){ printf("+"); fflush(stdout); }
997
998 idel = -idel ; /* reverse direction of ii */
999 }
1000 if( !be_quiet ){ printf("\n"); fflush(stdout); }
1001
1002 /*-- check diagonal for OK-ness --**/
1003
1004 atrace = 0.0 ;
1005 ii = 0 ;
1006 for( jj=0 ; jj < num_cols ; jj++ ){
1007 if( cov_mat[jj+jj*num_cols] <= 0.0 ){
1008 fprintf(stderr,"*** covariance diagonal (%d,%d) = %g\n",
1009 jj+1,jj+1,cov_mat[jj+jj*num_cols]) ;
1010 ii++ ;
1011 }
1012 atrace += cov_mat[jj+jj*num_cols] ;
1013 }
1014 if( ii > 0 ){ fprintf(stderr, "*** Warning %d zero or negative covariance on diagonals!\n", ii); }
1015
1016 if( !be_quiet ){ printf("--- covariance trace = %g\n",atrace); fflush(stdout); }
1017
1018 return(atrace);
1019 }
1020
1021 /*!
1022 Principal Component calculation by doing SVD on the
1023 covariance matrix of the data.
1024
1025 Based on code in 3dpc
1026
1027 data_mat is a matrix of M (num_cols) column vectors, each N (num_rows) elements
1028 long, and stored in a column major order.
1029 row_mask is a byte mask vector for data values to consider in data_mat
1030 If row_mask[i] then row i is considered in the calculations
1031 If NULL then all rows are considered.
1032 num_rows is the length N of each column in data_mat
1033 num_cols is the number of columns (M) (second dimension of data_mat)
1034 be_quiet
1035
1036 (HAVE YET TO MAKE THIS FUNCTION RETURN VALUES a la pca_fast3. I'll do it when
1037 I'll need it ZSS)
1038
1039 */
pca(float * data_mat,unsigned char * row_mask,int num_rows,int num_cols,int be_quiet)1040 void pca (float *data_mat, unsigned char * row_mask, int num_rows, int num_cols, int be_quiet) {
1041 double *aa=NULL, atrace, *wout, sum, *perc;
1042 int i, jj, ll, ii;
1043
1044 /* calculate covariance matrix */
1045 aa = (double *)malloc(sizeof(double)*num_cols*num_cols);
1046 wout = (double*)malloc(sizeof(double)*num_cols);
1047 atrace = covariance(data_mat, aa, row_mask, num_rows, num_cols, 0, 1, be_quiet);
1048
1049 /* calculate eigenvalues and vectors to be held in aa */
1050 symeig_double( num_cols , aa , wout ) ;
1051
1052 /* print results for now */
1053 sum = 0.0 ;
1054 perc = (double *) malloc( sizeof(double) * num_cols) ;
1055
1056 fprintf(stderr, "deal: Num. --Eigenvalue-- -Var.Fraction- -Cumul.Fract.-\n") ;
1057 for( jj=0 ; jj < num_cols ; jj++ ){
1058 ll = num_cols - 1-jj ; /* reversed order of eigensolution! */
1059 perc[jj] = wout[ll]/atrace ;
1060 sum += perc[jj] ;
1061 fprintf(stderr, "%4d %14.7g %14.7g %14.7g\n",
1062 jj+1 , wout[ll] , perc[jj] , sum ) ;
1063 }
1064
1065 /* now write the vectors */
1066 for (ii=0; ii< num_cols; ++ii) { /* for each row */
1067 for( jj=0 ; jj < num_cols ; jj++ ){ /* for each column */
1068 ll = num_cols - 1- jj ; /* reversed order of eigensolution! */
1069 fprintf(stderr, "%3.4f ",
1070 aa[ii+ll*num_cols] ) ;
1071 }
1072 fprintf(stderr, "\n"); fflush(stdout);
1073 }
1074
1075 free(perc); free(aa); free(wout);
1076 return;
1077 }
1078
1079 /*!
1080 Principal Component calculation by doing SVD on the
1081 covariance matrix of the data.
1082
1083 Optimized for matrices with 3 columns (x y z, typically)
1084 Based on code in 3dpc
1085
1086 data_mat is a matrix of 3 column vectors, each N (num_rows) elements
1087 long, and stored in a column major order.
1088 x0, x1, x2, x3, ... , y0, y1, ... , z0, z1, ... ,zN-1
1089
1090 num_rows is the length N of each column in data_mat
1091 be_quiet
1092 pca_vec is a matrix of 3 column vectors corresponding to the
1093 eigen values in pca_eig. NOTE: Storage is still column
1094 major order, like data_mat
1095 pca_eig is a vector of the 3 eigen values (sorted large to small)
1096
1097 Function returns the trace of the covariance matrix
1098 */
1099
pca_fast3(float * data_mat,int num_rows,int be_quiet,double * pca_mat,double * pca_eig)1100 double pca_fast3 (float *data_mat, int num_rows, int be_quiet, double *pca_mat, double *pca_eig) {
1101 double aa[9], atrace, wout[9], sum, perc[9];
1102 int i, jj, ll, ii, cnt;
1103
1104 #if DBG_PC
1105 fprintf(stderr, "data_mat=[ %f %f %f\n%f %f %f\n%f %f %f]\n",
1106 data_mat[0], data_mat[3], data_mat[6],
1107 data_mat[1], data_mat[4], data_mat[7],
1108 data_mat[2], data_mat[5], data_mat[8]);
1109 #endif
1110 /* calculate covariance */
1111 atrace = covariance(data_mat, aa, NULL, num_rows, 3, 0, 1, be_quiet);
1112
1113 /* calc eigs */
1114 symeig_3( aa , wout, 1 ) ;
1115
1116 /* now write the vectors */
1117 /* print results for now */
1118 sum = 0.0 ;
1119 #if DBG_PC
1120 fprintf(stderr, "deal3: Num. --Eigenvalue-- -Var.Fraction- -Cumul.Fract.-\n") ;
1121 #endif
1122 for( jj=0 ; jj < 3 ; jj++ ){
1123 ll = 3 - 1-jj ; /* reversed order of eigensolution! */
1124 pca_eig[jj] = wout[ll];
1125 #if DBG_PC
1126 perc[jj] = wout[ll]/atrace ;
1127 sum += perc[jj] ;
1128 fprintf(stderr, "%4d %14.7g %14.7g %14.7g\n",
1129 jj+1 , wout[ll] , perc[jj] , sum ) ;
1130 #endif
1131 }
1132
1133 cnt = 0;
1134 for (ii=0; ii< 3; ++ii) { /* for each row */
1135 for( jj=2 ; jj > -1 ; jj-- ){ /* for each column (reversed order) */
1136 ll = ii+jj*3;
1137 #if DBG_PC
1138 fprintf(stderr, "%3.4f ",
1139 aa[ll] ) ;
1140 #endif
1141 pca_mat[cnt] = aa[ll];
1142 ++cnt;
1143 }
1144 #if DBG_PC
1145 fprintf(stderr, "\n");
1146 #endif
1147 }
1148
1149 return(atrace);
1150
1151 }
1152
1153 /*---------------------------------------------------------------------------*/
1154
1155 #undef A
1156 #undef U
1157 #undef V
1158 #undef pseps
1159
1160 #define A(i,j) aa[(i)+(j)*m]
1161 #define U(i,j) uu[(i)+(j)*m]
1162 #define V(i,j) vv[(i)+(j)*n]
1163 #define pseps 5.e-7
1164
1165 #undef DSDBUG
1166
svd_desingularize(int m,int n,double * aa)1167 int svd_desingularize( int m , int n , double *aa )
1168 {
1169 double *uu , *vv , *ww , sum , smax ;
1170 int i,j,k , nfix ;
1171
1172 if( aa == NULL || m < 1 || n < 1 ) return -1 ;
1173
1174 ww = (double *)malloc(sizeof(double)*n) ;
1175 uu = (double *)malloc(sizeof(double)*m*n) ;
1176 vv = (double *)malloc(sizeof(double)*n*n) ;
1177
1178 svd_double( m , n , aa , ww , uu , vv ) ;
1179
1180 /* fix singular values */
1181
1182 smax = ww[0] ;
1183 for( i=1 ; i < n ; i++ ) if( ww[i] > smax ) smax = ww[i] ;
1184 if( smax < 0.0 ){ /* should not happen */
1185 free(vv) ; free(uu) ; free(ww) ; return -1 ;
1186 } else if( smax == 0.0 ){ /* only if input is all zero */
1187 smax = 1.0 ;
1188 }
1189
1190 #ifdef DSDBUG
1191 fprintf(stderr,"============= Desingularization ==============\n") ;
1192 fprintf(stderr,"ww before =") ;
1193 for( i=0 ; i < n ; i++ ) fprintf(stderr," %g",ww[i]) ;
1194 fprintf(stderr,"\n") ;
1195 #endif
1196
1197 smax = pseps * smax ;
1198 for( nfix=i=0 ; i < n ; i++ ){
1199 if( ww[i] < 0.0 ){ ww[i] = smax ; nfix++; }
1200 else if( ww[i] < 2.0*smax ){ ww[i] = smax+0.25*ww[i]*ww[i]/smax; nfix++; }
1201 }
1202
1203 #ifdef DSDBUG
1204 fprintf(stderr,"ww after =") ;
1205 for( i=0 ; i < n ; i++ ) fprintf(stderr," %g",ww[i]) ;
1206 fprintf(stderr,"\n") ;
1207 #endif
1208
1209 if( nfix == 0 ){ /* no fix needed */
1210 #ifdef DSDBUG
1211 fprintf(stderr,".......... no desingularization needed ..........\n") ;
1212 #endif
1213 free(vv) ; free(uu) ; free(ww) ; return 0 ;
1214 }
1215
1216 #ifdef DSDBUG
1217 fprintf(stderr,"U\n") ;
1218 for( i=0 ; i < m ; i++ ){
1219 for( j=0 ; j < n ; j++ ) fprintf(stderr," %g",U(i,j)) ;
1220 fprintf(stderr,"\n") ;
1221 }
1222 fprintf(stderr,"V\n") ;
1223 for( i=0 ; i < n ; i++ ){
1224 for( j=0 ; j < n ; j++ ) fprintf(stderr," %g",V(i,j)) ;
1225 fprintf(stderr,"\n") ;
1226 }
1227 #endif
1228
1229 /* otherwise, recompute [a] matrix from fixed SVD */
1230
1231 for( j=0 ; j < n ; j++ ){
1232 for( i=0 ; i < m ; i++ ){
1233 sum = 0.0 ;
1234 for( k=0 ; k < n ; k++ ) sum += U(i,k)*V(j,k)*ww[k] ;
1235 #ifdef DSDBUG
1236 fprintf(stderr,"A(%d,%d): %g -> %g\n",i,j,A(i,j),sum) ;
1237 #endif
1238 A(i,j) = sum ;
1239 }}
1240 #ifdef DSDBUG
1241 fprintf(stderr,"............ end desingularization ..............\n") ;
1242 #endif
1243
1244 free(vv) ; free(uu) ; free(ww) ; return nfix ;
1245 }
1246