1 #ifndef _I_ALREADY_HAVE_SVDLIB_
2 #define _I_ALREADY_HAVE_SVDLIB_
3 /*
4 Copyright © 2002, University of Tennessee Research Foundation.
5 All rights reserved.
6
7 Redistribution and use in source and binary forms, with or without
8 modification, are permitted provided that the following conditions are met:
9
10 * Redistributions of source code must retain the above copyright notice, this
11 list of conditions and the following disclaimer.
12
13 Redistributions in binary form must reproduce the above copyright notice,
14 this list of conditions and the following disclaimer in the documentation
15 and/or other materials provided with the distribution.
16
17 * Neither the name of the University of Tennessee nor the names of its
18 contributors may be used to endorse or promote products derived from this
19 software without specific prior written permission.
20
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
25 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 POSSIBILITY OF SUCH DAMAGE.
32 */
33
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <stdarg.h>
37 #include <string.h>
38 #include <errno.h>
39 #include <math.h>
40 #include <sys/types.h>
41 #include <sys/stat.h>
42 #include <netinet/in.h>
43 #include <fcntl.h>
44 #include "svdlib.h"
45
46 #define BUNZIP2 "bzip2 -d"
47 #define BZIP2 "bzip2 -1"
48 #define UNZIP "gzip -d"
49 #define ZIP "gzip -1"
50 #define COMPRESS "compress"
51
52 #define MAX_FILENAME 512
53 #define MAX_PIPES 64
54 static FILE *Pipe[MAX_PIPES];
55 static int numPipes = 0;
56
svd_longArray(long size,char empty,char * name)57 long *svd_longArray(long size, char empty, char *name) {
58 long *a;
59 if (empty) a = (long *) calloc(size, sizeof(long));
60 else a = (long *) malloc(size * sizeof(long));
61 if (!a) {
62 perror(name);
63 /* exit(errno); */
64 }
65 return a;
66 }
67
svd_doubleArray(long size,char empty,char * name)68 double *svd_doubleArray(long size, char empty, char *name) {
69 double *a;
70 if (empty) a = (double *) calloc(size, sizeof(double));
71 else a = (double *) malloc(size * sizeof(double));
72 if (!a) {
73 perror(name);
74 /* exit(errno); */
75 }
76 return a;
77 }
78
svd_beep(void)79 void svd_beep(void) {
80 fputc('\a', stderr);
81 fflush(stderr);
82 }
83
svd_debug(char * fmt,...)84 void svd_debug(char *fmt, ...) {
85 va_list args;
86 va_start(args, fmt);
87 vfprintf(stderr, fmt, args);
88 va_end(args);
89 }
90
svd_error(char * fmt,...)91 void svd_error(char *fmt, ...) {
92 va_list args;
93 va_start(args, fmt);
94 svd_beep();
95 fprintf(stderr, "ERROR: ");
96 vfprintf(stderr, fmt, args);
97 fprintf(stderr, "\n");
98 va_end(args);
99 }
100
svd_fatalError(char * fmt,...)101 void svd_fatalError(char *fmt, ...) {
102 va_list args;
103 va_start(args, fmt);
104 svd_beep();
105 fprintf(stderr, "ERROR: ");
106 vfprintf(stderr, fmt, args);
107 fprintf(stderr, "\a\n");
108 va_end(args);
109 exit(1);
110 }
111
registerPipe(FILE * p)112 static void registerPipe(FILE *p) {
113 if (numPipes >= MAX_PIPES) svd_error("Too many pipes open");
114 Pipe[numPipes++] = p;
115 }
116
isPipe(FILE * p)117 static char isPipe(FILE *p) {
118 int i;
119 for (i = 0; i < numPipes && Pipe[i] != p; i++);
120 if (i == numPipes) return FALSE;
121 Pipe[i] = Pipe[--numPipes];
122 return TRUE;
123 }
124
openPipe(char * pipeName,char * mode)125 static FILE *openPipe(char *pipeName, char *mode) {
126 FILE *pipe;
127 fflush(stdout);
128 if ((pipe = popen(pipeName, mode))) registerPipe(pipe);
129 return pipe;
130 }
131
readZippedFile(char * command,char * fileName)132 static FILE *readZippedFile(char *command, char *fileName) {
133 char buf[MAX_FILENAME];
134 sprintf(buf, "%s < %s 2>/dev/null", command, fileName);
135 return openPipe(buf, "r");
136 }
137
svd_fatalReadFile(char * filename)138 FILE *svd_fatalReadFile(char *filename) {
139 FILE *file;
140 if (!(file = svd_readFile(filename)))
141 svd_fatalError("couldn't read the file %s", filename);
142 return file;
143 }
144
stringEndsIn(char * s,char * t)145 static int stringEndsIn(char *s, char *t) {
146 int ls = strlen(s);
147 int lt = strlen(t);
148 if (ls < lt) return FALSE;
149 return (strcmp(s + ls - lt, t)) ? FALSE : TRUE;
150 }
151
152 /* Will silently return NULL if file couldn't be opened */
svd_readFile(char * fileName)153 FILE *svd_readFile(char *fileName) {
154 char fileBuf[MAX_FILENAME];
155 struct stat statbuf;
156
157 /* Special file name */
158 if (!strcmp(fileName, "-"))
159 return stdin;
160
161 /* If it is a pipe */
162 if (fileName[0] == '|')
163 return openPipe(fileName + 1, "r");
164
165 /* Check if already ends in .gz or .Z and assume compressed */
166 if (stringEndsIn(fileName, ".gz") || stringEndsIn(fileName, ".Z")) {
167 if (!stat(fileName, &statbuf))
168 return readZippedFile(UNZIP, fileName);
169 return NULL;
170 }
171 /* Check if already ends in .bz or .bz2 and assume compressed */
172 if (stringEndsIn(fileName, ".bz") || stringEndsIn(fileName, ".bz2")) {
173 if (!stat(fileName, &statbuf))
174 return readZippedFile(BUNZIP2, fileName);
175 return NULL;
176 }
177 /* Try just opening normally */
178 if (!stat(fileName, &statbuf))
179 return fopen(fileName, "r");
180 /* Try adding .gz */
181 sprintf(fileBuf, "%s.gz", fileName);
182 if (!stat(fileBuf, &statbuf))
183 return readZippedFile(UNZIP, fileBuf);
184 /* Try adding .Z */
185 sprintf(fileBuf, "%s.Z", fileName);
186 if (!stat(fileBuf, &statbuf))
187 return readZippedFile(UNZIP, fileBuf);
188 /* Try adding .bz2 */
189 sprintf(fileBuf, "%s.bz2", fileName);
190 if (!stat(fileBuf, &statbuf))
191 return readZippedFile(BUNZIP2, fileBuf);
192 /* Try adding .bz */
193 sprintf(fileBuf, "%s.bz", fileName);
194 if (!stat(fileBuf, &statbuf))
195 return readZippedFile(BUNZIP2, fileBuf);
196
197 return NULL;
198 }
199
writeZippedFile(char * fileName,char append)200 static FILE *writeZippedFile(char *fileName, char append) {
201 char buf[MAX_FILENAME];
202 const char *op = (append) ? ">>" : ">";
203 if (stringEndsIn(fileName, ".bz2") || stringEndsIn(fileName, ".bz"))
204 sprintf(buf, "%s %s \"%s\"", BZIP2, op, fileName);
205 else if (stringEndsIn(fileName, ".Z"))
206 sprintf(buf, "%s %s \"%s\"", COMPRESS, op, fileName);
207 else
208 sprintf(buf, "%s %s \"%s\"", ZIP, op, fileName);
209 return openPipe(buf, "w");
210 }
211
svd_writeFile(char * fileName,char append)212 FILE *svd_writeFile(char *fileName, char append) {
213 /* Special file name */
214 if (!strcmp(fileName, "-"))
215 return stdout;
216
217 /* If it is a pipe */
218 if (fileName[0] == '|')
219 return openPipe(fileName + 1, "w");
220
221 /* Check if ends in .gz, .Z, .bz, .bz2 */
222 if (stringEndsIn(fileName, ".gz") || stringEndsIn(fileName, ".Z") ||
223 stringEndsIn(fileName, ".bz") || stringEndsIn(fileName, ".bz2"))
224 return writeZippedFile(fileName, append);
225 return (append) ? fopen(fileName, "a") : fopen(fileName, "w");
226 }
227
228 /* Could be a file or a stream. */
svd_closeFile(FILE * file)229 void svd_closeFile(FILE *file) {
230 if (file == stdin || file == stdout) return;
231 if (isPipe(file)) pclose(file);
232 else fclose(file);
233 }
234
235
svd_readBinInt(FILE * file,int * val)236 char svd_readBinInt(FILE *file, int *val) {
237 int x;
238 if (fread(&x, sizeof(int), 1, file) == 1) {
239 *val = ntohl(x);
240 return FALSE;
241 }
242 return TRUE;
243 }
244
245 /* This reads a float in network order and converts to a real in host order. */
svd_readBinFloat(FILE * file,float * val)246 char svd_readBinFloat(FILE *file, float *val) {
247 int x;
248 float y;
249 if (fread(&x, sizeof(int), 1, file) == 1) {
250 x = ntohl(x);
251 y = *((float *) &x);
252 *val = y;
253 return FALSE;
254 }
255 return TRUE;
256 }
257
svd_writeBinInt(FILE * file,int x)258 char svd_writeBinInt(FILE *file, int x) {
259 int y = htonl(x);
260 if (fwrite(&y, sizeof(int), 1, file) != 1) return TRUE;
261 return FALSE;
262 }
263
264 /* This takes a real in host order and writes a float in network order. */
svd_writeBinFloat(FILE * file,float r)265 char svd_writeBinFloat(FILE *file, float r) {
266 int y = htonl(*((int *) &r));
267 if (fwrite(&y, sizeof(int), 1, file) != 1) return TRUE;
268 return FALSE;
269 }
270
271
272 /**************************************************************
273 * returns |a| if b is positive; else fsign returns -|a| *
274 **************************************************************/
svd_fsign(double a,double b)275 double svd_fsign(double a, double b) {
276 if ((a>=0.0 && b>=0.0) || (a<0.0 && b<0.0))return(a);
277 else return -a;
278 }
279
280 /**************************************************************
281 * returns the larger of two double precision numbers *
282 **************************************************************/
svd_dmax(double a,double b)283 double svd_dmax(double a, double b) {
284 return (a > b) ? a : b;
285 }
286
287 /**************************************************************
288 * returns the smaller of two double precision numbers *
289 **************************************************************/
svd_dmin(double a,double b)290 double svd_dmin(double a, double b) {
291 return (a < b) ? a : b;
292 }
293
294 /**************************************************************
295 * returns the larger of two integers *
296 **************************************************************/
svd_imax(long a,long b)297 long svd_imax(long a, long b) {
298 return (a > b) ? a : b;
299 }
300
301 /**************************************************************
302 * returns the smaller of two integers *
303 **************************************************************/
svd_imin(long a,long b)304 long svd_imin(long a, long b) {
305 return (a < b) ? a : b;
306 }
307
308 /**************************************************************
309 * Function scales a vector by a constant. *
310 * Based on Fortran-77 routine from Linpack by J. Dongarra *
311 **************************************************************/
svd_dscal(long n,double da,double * dx,long incx)312 void svd_dscal(long n, double da, double *dx, long incx) {
313 long i;
314
315 if (n <= 0 || incx == 0) return;
316 if (incx < 0) dx += (-n+1) * incx;
317 for (i=0; i < n; i++) {
318 *dx *= da;
319 dx += incx;
320 }
321 return;
322 }
323
324 /**************************************************************
325 * function scales a vector by a constant. *
326 * Based on Fortran-77 routine from Linpack by J. Dongarra *
327 **************************************************************/
svd_datx(long n,double da,double * dx,long incx,double * dy,long incy)328 void svd_datx(long n, double da, double *dx, long incx, double *dy, long incy) {
329 long i;
330
331 if (n <= 0 || incx == 0 || incy == 0 || da == 0.0) return;
332 if (incx == 1 && incy == 1)
333 for (i=0; i < n; i++) *dy++ = da * (*dx++);
334
335 else {
336 if (incx < 0) dx += (-n+1) * incx;
337 if (incy < 0) dy += (-n+1) * incy;
338 for (i=0; i < n; i++) {
339 *dy = da * (*dx);
340 dx += incx;
341 dy += incy;
342 }
343 }
344 return;
345 }
346
347 /**************************************************************
348 * Function copies a vector x to a vector y *
349 * Based on Fortran-77 routine from Linpack by J. Dongarra *
350 **************************************************************/
svd_dcopy(long n,double * dx,long incx,double * dy,long incy)351 void svd_dcopy(long n, double *dx, long incx, double *dy, long incy) {
352 long i;
353
354 if (n <= 0 || incx == 0 || incy == 0) return;
355 if (incx == 1 && incy == 1)
356 for (i=0; i < n; i++) *dy++ = *dx++;
357
358 else {
359 if (incx < 0) dx += (-n+1) * incx;
360 if (incy < 0) dy += (-n+1) * incy;
361 for (i=0; i < n; i++) {
362 *dy = *dx;
363 dx += incx;
364 dy += incy;
365 }
366 }
367 return;
368 }
369
370 /**************************************************************
371 * Function forms the dot product of two vectors. *
372 * Based on Fortran-77 routine from Linpack by J. Dongarra *
373 **************************************************************/
svd_ddot(long n,double * dx,long incx,double * dy,long incy)374 double svd_ddot(long n, double *dx, long incx, double *dy, long incy) {
375 long i;
376 double dot_product;
377
378 if (n <= 0 || incx == 0 || incy == 0) return(0.0);
379 dot_product = 0.0;
380 if (incx == 1 && incy == 1)
381 for (i=0; i < n; i++) dot_product += (*dx++) * (*dy++);
382 else {
383 if (incx < 0) dx += (-n+1) * incx;
384 if (incy < 0) dy += (-n+1) * incy;
385 for (i=0; i < n; i++) {
386 dot_product += (*dx) * (*dy);
387 dx += incx;
388 dy += incy;
389 }
390 }
391 return(dot_product);
392 }
393
394 /**************************************************************
395 * Constant times a vector plus a vector *
396 * Based on Fortran-77 routine from Linpack by J. Dongarra *
397 **************************************************************/
svd_daxpy(long n,double da,double * dx,long incx,double * dy,long incy)398 void svd_daxpy (long n, double da, double *dx, long incx, double *dy, long incy) {
399 long i;
400
401 if (n <= 0 || incx == 0 || incy == 0 || da == 0.0) return;
402 if (incx == 1 && incy == 1)
403 for (i=0; i < n; i++) {
404 *dy += da * (*dx++);
405 dy++;
406 }
407 else {
408 if (incx < 0) dx += (-n+1) * incx;
409 if (incy < 0) dy += (-n+1) * incy;
410 for (i=0; i < n; i++) {
411 *dy += da * (*dx);
412 dx += incx;
413 dy += incy;
414 }
415 }
416 return;
417 }
418
419 /*********************************************************************
420 * Function sorts array1 and array2 into increasing order for array1 *
421 *********************************************************************/
svd_dsort2(long igap,long n,double * array1,double * array2)422 void svd_dsort2(long igap, long n, double *array1, double *array2) {
423 double temp;
424 long i, j, index;
425
426 if (!igap) return;
427 else {
428 for (i = igap; i < n; i++) {
429 j = i - igap;
430 index = i;
431 while (j >= 0 && array1[j] > array1[index]) {
432 temp = array1[j];
433 array1[j] = array1[index];
434 array1[index] = temp;
435 temp = array2[j];
436 array2[j] = array2[index];
437 array2[index] = temp;
438 j -= igap;
439 index = j + igap;
440 }
441 }
442 }
443 svd_dsort2(igap/2,n,array1,array2);
444 }
445
446 /**************************************************************
447 * Function interchanges two vectors *
448 * Based on Fortran-77 routine from Linpack by J. Dongarra *
449 **************************************************************/
svd_dswap(long n,double * dx,long incx,double * dy,long incy)450 void svd_dswap(long n, double *dx, long incx, double *dy, long incy) {
451 long i;
452 double dtemp;
453
454 if (n <= 0 || incx == 0 || incy == 0) return;
455 if (incx == 1 && incy == 1) {
456 for (i=0; i < n; i++) {
457 dtemp = *dy;
458 *dy++ = *dx;
459 *dx++ = dtemp;
460 }
461 }
462 else {
463 if (incx < 0) dx += (-n+1) * incx;
464 if (incy < 0) dy += (-n+1) * incy;
465 for (i=0; i < n; i++) {
466 dtemp = *dy;
467 *dy = *dx;
468 *dx = dtemp;
469 dx += incx;
470 dy += incy;
471 }
472 }
473 }
474
475 /*****************************************************************
476 * Function finds the index of element having max. absolute value*
477 * based on FORTRAN 77 routine from Linpack by J. Dongarra *
478 *****************************************************************/
svd_idamax(long n,double * dx,long incx)479 long svd_idamax(long n, double *dx, long incx) {
480 long ix,i,imax;
481 double dtemp, dmax;
482
483 if (n < 1) return(-1);
484 if (n == 1) return(0);
485 if (incx == 0) return(-1);
486
487 if (incx < 0) ix = (-n+1) * incx;
488 else ix = 0;
489 imax = ix;
490 dx += ix;
491 dmax = fabs(*dx);
492 for (i=1; i < n; i++) {
493 ix += incx;
494 dx += incx;
495 dtemp = fabs(*dx);
496 if (dtemp > dmax) {
497 dmax = dtemp;
498 imax = ix;
499 }
500 }
501 return(imax);
502 }
503
504 /**************************************************************
505 * multiplication of matrix B by vector x, where B = A'A, *
506 * and A is nrow by ncol (nrow >> ncol). Hence, B is of order *
507 * n = ncol (y stores product vector). *
508 **************************************************************/
svd_opb(SMat A,double * x,double * y,double * temp)509 void svd_opb(SMat A, double *x, double *y, double *temp) {
510 long i, j, end;
511 long *pointr = A->pointr, *rowind = A->rowind;
512 double *value = A->value;
513 long n = A->cols;
514
515 #ifndef USE_OMP
516 SVDCount[SVD_MXV] += 2;
517 #endif
518 memset(y, 0, n * sizeof(double));
519 for (i = 0; i < A->rows; i++) temp[i] = 0.0;
520
521 for (i = 0; i < A->cols; i++) {
522 end = pointr[i+1];
523 for (j = pointr[i]; j < end; j++)
524 temp[rowind[j]] += value[j] * (*x);
525 x++;
526 }
527 for (i = 0; i < A->cols; i++) {
528 end = pointr[i+1];
529 for (j = pointr[i]; j < end; j++)
530 *y += value[j] * temp[rowind[j]];
531 y++;
532 }
533 return;
534 }
535
536 /***********************************************************
537 * multiplication of matrix A by vector x, where A is *
538 * nrow by ncol (nrow >> ncol). y stores product vector. *
539 ***********************************************************/
svd_opa(SMat A,double * x,double * y)540 void svd_opa(SMat A, double *x, double *y) {
541 long end, i, j;
542 long *pointr = A->pointr, *rowind = A->rowind;
543 double *value = A->value;
544
545 #ifndef USE_OMP
546 SVDCount[SVD_MXV]++;
547 #endif
548 memset(y, 0, A->rows * sizeof(double));
549
550 for (i = 0; i < A->cols; i++) {
551 end = pointr[i+1];
552 for (j = pointr[i]; j < end; j++)
553 y[rowind[j]] += value[j] * x[i];
554 }
555 return;
556 }
557
558
559 /***********************************************************************
560 * *
561 * random() *
562 * (double precision) *
563 ***********************************************************************/
564 /***********************************************************************
565
566 Description
567 -----------
568
569 This is a translation of a Fortran-77 uniform random number
570 generator. The code is based on theory and suggestions given in
571 D. E. Knuth (1969), vol 2. The argument to the function should
572 be initialized to an arbitrary integer prior to the first call to
573 random. The calling program should not alter the value of the
574 argument between subsequent calls to random. Random returns values
575 within the interval (0,1).
576
577
578 Arguments
579 ---------
580
581 (input)
582 iy an integer seed whose value must not be altered by the caller
583 between subsequent calls
584
585 (output)
586 random a double precision random number between (0,1)
587
588 ***********************************************************************/
svd_random2(long * iy)589 double svd_random2(long *iy) {
590 static long m2 = 0;
591 static long ia, ic, mic;
592 static double halfm, s;
593
594 /* If first entry, compute (max int) / 2 */
595 { if (!m2) {
596 m2 = 1 << (8 * (int)sizeof(int) - 2);
597 halfm = m2;
598
599 /* compute multiplier and increment for linear congruential method */
600 ia = 8 * (long)(halfm * atan(1.0) / 8.0) + 5;
601 ic = 2 * (long)(halfm * (0.5 - sqrt(3.0)/6.0)) + 1;
602 mic = (m2-ic) + m2;
603
604 /* s is the scale factor for converting to floating point */
605 s = 0.5 / halfm;
606 }}
607 if( iy == NULL ) return 0.0 ; /* RWCox */
608
609 /* compute next random number */
610 *iy = *iy * ia;
611
612 /* for computers which do not allow integer overflow on addition */
613 if (*iy > mic) *iy = (*iy - m2) - m2;
614
615 *iy = *iy + ic;
616
617 /* for computers whose word length for addition is greater than
618 * for multiplication */
619 if (*iy / 2 > m2) *iy = (*iy - m2) - m2;
620
621 /* for computers whose integer overflow affects the sign bit */
622 if (*iy < 0) *iy = (*iy + m2) + m2;
623
624 return((double)(*iy) * s);
625 }
626
627 /**************************************************************
628 * *
629 * Function finds sqrt(a^2 + b^2) without overflow or *
630 * destructive underflow. *
631 * *
632 **************************************************************/
633 /**************************************************************
634
635 Funtions used
636 -------------
637
638 UTILITY dmax, dmin
639
640 **************************************************************/
svd_pythag(double a,double b)641 double svd_pythag(double a, double b) {
642 double p, r, s, t, u, temp;
643
644 p = svd_dmax(fabs(a), fabs(b));
645 if (p != 0.0) {
646 temp = svd_dmin(fabs(a), fabs(b)) / p;
647 r = temp * temp;
648 t = 4.0 + r;
649 while (t != 4.0) {
650 s = r / t;
651 u = 1.0 + 2.0 * s;
652 p *= u;
653 temp = s / u;
654 r *= temp * temp;
655 t = 4.0 + r;
656 }
657 }
658 return(p);
659 }
660
661 /*
662 Copyright © 2002, University of Tennessee Research Foundation.
663 All rights reserved.
664
665 Redistribution and use in source and binary forms, with or without
666 modification, are permitted provided that the following conditions are met:
667
668 * Redistributions of source code must retain the above copyright notice, this
669 list of conditions and the following disclaimer.
670
671 Redistributions in binary form must reproduce the above copyright notice,
672 this list of conditions and the following disclaimer in the documentation
673 and/or other materials provided with the distribution.
674
675 * Neither the name of the University of Tennessee nor the names of its
676 contributors may be used to endorse or promote products derived from this
677 software without specific prior written permission.
678
679 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
680 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
681 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
682 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
683 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
684 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
685 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
686 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
687 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
688 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
689 POSSIBILITY OF SUCH DAMAGE.
690 */
691
692 char *SVDVersion = "1.4";
693 long SVDVerbosity = 0;
694 long SVDCount[SVD_COUNTERS];
695
696 #ifndef USE_OMP
svdResetCounters(void)697 void svdResetCounters(void) {
698 int i;
699 for (i = 0; i < SVD_COUNTERS; i++)
700 SVDCount[i] = 0;
701 }
702 #else
svdResetCounters(void)703 void svdResetCounters(void){ return; }
704 #endif
705
706 /********************************* Allocation ********************************/
707
708 /* Row major order. Rows are vectors that are consecutive in memory. Matrix
709 is initialized to empty. */
svdNewDMat(int rows,int cols)710 DMat svdNewDMat(int rows, int cols) {
711 int i;
712 DMat D = (DMat) malloc(sizeof(struct dmat));
713 if (!D) {perror("svdNewDMat"); return NULL;}
714 D->rows = rows;
715 D->cols = cols;
716
717 D->value = (double **) malloc(rows * sizeof(double *));
718 if (!D->value) {SAFE_FREE(D); return NULL;}
719
720 D->value[0] = (double *) calloc(rows * cols, sizeof(double));
721 if (!D->value[0]) {SAFE_FREE(D->value); SAFE_FREE(D); return NULL;}
722
723 for (i = 1; i < rows; i++) D->value[i] = D->value[i-1] + cols;
724 return D;
725 }
726
svdFreeDMat(DMat D)727 void svdFreeDMat(DMat D) {
728 if (!D) return;
729 SAFE_FREE(D->value[0]);
730 SAFE_FREE(D->value);
731 free(D);
732 }
733
734
svdNewSMat(int rows,int cols,int vals)735 SMat svdNewSMat(int rows, int cols, int vals) {
736 SMat S = (SMat) calloc(1, sizeof(struct smat));
737 if (!S) {perror("svdNewSMat"); return NULL;}
738 S->rows = rows;
739 S->cols = cols;
740 S->vals = vals;
741 S->pointr = svd_longArray(cols + 1, TRUE, "svdNewSMat: pointr");
742 if (!S->pointr) {svdFreeSMat(S); return NULL;}
743 S->rowind = svd_longArray(vals, FALSE, "svdNewSMat: rowind");
744 if (!S->rowind) {svdFreeSMat(S); return NULL;}
745 S->value = svd_doubleArray(vals, FALSE, "svdNewSMat: value");
746 if (!S->value) {svdFreeSMat(S); return NULL;}
747 return S;
748 }
749
svdFreeSMat(SMat S)750 void svdFreeSMat(SMat S) {
751 if (!S) return;
752 SAFE_FREE(S->pointr);
753 SAFE_FREE(S->rowind);
754 SAFE_FREE(S->value);
755 free(S);
756 }
757
758
759 /* Creates an empty SVD record */
svdNewSVDRec(void)760 SVDRec svdNewSVDRec(void) {
761 SVDRec R = (SVDRec) calloc(1, sizeof(struct svdrec));
762 if (!R) {perror("svdNewSVDRec"); return NULL;}
763 return R;
764 }
765
766 /* Frees an svd rec and all its contents. */
svdFreeSVDRec(SVDRec R)767 void svdFreeSVDRec(SVDRec R) {
768 if (!R) return;
769 if (R->Ut) svdFreeDMat(R->Ut);
770 if (R->S) SAFE_FREE(R->S);
771 if (R->Vt) svdFreeDMat(R->Vt);
772 free(R);
773 }
774
775
776 /**************************** Conversion *************************************/
777
778 /* Converts a sparse matrix to a dense one (without affecting the former) */
svdConvertStoD(SMat S)779 DMat svdConvertStoD(SMat S) {
780 int i, c;
781 DMat D = svdNewDMat(S->rows, S->cols);
782 if (!D) {
783 svd_error("svdConvertStoD: failed to allocate D");
784 return NULL;
785 }
786 for (i = 0, c = 0; i < S->vals; i++) {
787 while (S->pointr[c + 1] <= i) c++;
788 D->value[S->rowind[i]][c] = S->value[i];
789 }
790 return D;
791 }
792
793 /* Converts a dense matrix to a sparse one (without affecting the dense one) */
svdConvertDtoS(DMat D)794 SMat svdConvertDtoS(DMat D) {
795 SMat S;
796 int i, j, n;
797 for (i = 0, n = 0; i < D->rows; i++)
798 for (j = 0; j < D->cols; j++)
799 if (D->value[i][j] != 0) n++;
800
801 S = svdNewSMat(D->rows, D->cols, n);
802 if (!S) {
803 svd_error("svdConvertDtoS: failed to allocate S");
804 return NULL;
805 }
806 for (j = 0, n = 0; j < D->cols; j++) {
807 S->pointr[j] = n;
808 for (i = 0; i < D->rows; i++)
809 if (D->value[i][j] != 0) {
810 S->rowind[n] = i;
811 S->value[n] = D->value[i][j];
812 n++;
813 }
814 }
815 S->pointr[S->cols] = S->vals;
816 return S;
817 }
818
819 /* Transposes a dense matrix. */
svdTransposeD(DMat D)820 DMat svdTransposeD(DMat D) {
821 int r, c;
822 DMat N = svdNewDMat(D->cols, D->rows);
823 for (r = 0; r < D->rows; r++)
824 for (c = 0; c < D->cols; c++)
825 N->value[c][r] = D->value[r][c];
826 return N;
827 }
828
829 /* Efficiently transposes a sparse matrix. */
svdTransposeS(SMat S)830 SMat svdTransposeS(SMat S) {
831 int r, c, i, j;
832 SMat N = svdNewSMat(S->cols, S->rows, S->vals);
833 /* Count number nz in each row. */
834 for (i = 0; i < S->vals; i++)
835 N->pointr[S->rowind[i]]++;
836 /* Fill each cell with the starting point of the previous row. */
837 N->pointr[S->rows] = S->vals - N->pointr[S->rows - 1];
838 for (r = S->rows - 1; r > 0; r--)
839 N->pointr[r] = N->pointr[r+1] - N->pointr[r-1];
840 N->pointr[0] = 0;
841 /* Assign the new columns and values. */
842 for (c = 0, i = 0; c < S->cols; c++) {
843 for (; i < S->pointr[c+1]; i++) {
844 r = S->rowind[i];
845 j = N->pointr[r+1]++;
846 N->rowind[j] = c;
847 N->value[j] = S->value[i];
848 }
849 }
850 return N;
851 }
852
853
854 /**************************** Input/Output ***********************************/
855
svdWriteDenseArray(double * a,int n,char * filename,char binary)856 void svdWriteDenseArray(double *a, int n, char *filename, char binary) {
857 int i;
858 FILE *file = svd_writeFile(filename, FALSE);
859 if (!file){
860 svd_error("svdWriteDenseArray: failed to write %s", filename);
861 return ;
862 }
863 if (binary) {
864 svd_writeBinInt(file, n);
865 for (i = 0; i < n; i++)
866 svd_writeBinFloat(file, (float) a[i]);
867 } else {
868 fprintf(file, "%d\n", n);
869 for (i = 0; i < n; i++)
870 fprintf(file, "%g ", a[i]);
871 fprintf(file,"\n") ;
872 }
873 svd_closeFile(file);
874 }
875
svdLoadDenseArray(char * filename,int * np,char binary)876 double *svdLoadDenseArray(char *filename, int *np, char binary) {
877 int i, n;
878 double *a;
879
880 FILE *file = svd_readFile(filename);
881 if (!file) {
882 svd_error("svdLoadDenseArray: failed to read %s", filename);
883 return NULL;
884 }
885 if (binary) {
886 svd_readBinInt(file, np);
887 } else if (fscanf(file, " %d", np) != 1) {
888 svd_error("svdLoadDenseArray: error reading %s", filename);
889 svd_closeFile(file);
890 return NULL;
891 }
892 n = *np;
893 a = svd_doubleArray(n, FALSE, "svdLoadDenseArray: a");
894 if (!a) return NULL;
895 if (binary) {
896 float f=0.0f;
897 for (i = 0; i < n; i++) {
898 svd_readBinFloat(file, &f);
899 a[i] = f;
900 }
901 } else {
902 for (i = 0; i < n; i++) {
903 if (fscanf(file, " %lf\n", a + i) != 1) {
904 svd_error("svdLoadDenseArray: error reading %s", filename);
905 break;
906 }
907 }
908 }
909 svd_closeFile(file);
910 return a;
911 }
912
913
914 /* File format has a funny header, then first entry index per column, then the
915 row for each entry, then the value for each entry. Indices count from 1.
916 Assumes A is initialized. */
svdLoadSparseTextHBFile(FILE * file)917 static SMat svdLoadSparseTextHBFile(FILE *file) {
918 char line[128];
919 long i, x, rows, cols, vals, num_mat;
920 SMat S;
921 /* Skip the header line: */
922 if (!fgets(line, 128, file));
923 /* Skip the line giving the number of lines in this file: */
924 if (!fgets(line, 128, file));
925 /* Read the line with useful dimensions: */
926 if (fscanf(file, "%*s%ld%ld%ld%ld\n",
927 &rows, &cols, &vals, &num_mat) != 4) {
928 svd_error("svdLoadSparseTextHBFile: bad file format on line 3");
929 return NULL;
930 }
931 if (num_mat != 0) {
932 svd_error("svdLoadSparseTextHBFile: I don't know how to handle a file "
933 "with elemental matrices (last entry on header line 3)");
934 return NULL;
935 }
936 /* Skip the line giving the formats: */
937 if (!fgets(line, 128, file));
938
939 S = svdNewSMat(rows, cols, vals);
940 if (!S) return NULL;
941
942 /* Read column pointers. */
943 for (i = 0; i <= S->cols; i++) {
944 if (fscanf(file, " %ld", &x) != 1) {
945 svd_error("svdLoadSparseTextHBFile: error reading pointr %d", i);
946 return NULL;
947 }
948 S->pointr[i] = x - 1;
949 }
950 S->pointr[S->cols] = S->vals;
951
952 /* Read row indices. */
953 for (i = 0; i < S->vals; i++) {
954 if (fscanf(file, " %ld", &x) != 1) {
955 svd_error("svdLoadSparseTextHBFile: error reading rowind %d", i);
956 return NULL;
957 }
958 S->rowind[i] = x - 1;
959 }
960 for (i = 0; i < S->vals; i++)
961 if (fscanf(file, " %lf", S->value + i) != 1) {
962 svd_error("svdLoadSparseTextHBFile: error reading value %d", i);
963 return NULL;
964 }
965 return S;
966 }
967
svdWriteSparseTextHBFile(SMat S,FILE * file)968 static void svdWriteSparseTextHBFile(SMat S, FILE *file) {
969 int i;
970 long col_lines = ((S->cols + 1) / 8) + (((S->cols + 1) % 8) ? 1 : 0);
971 long row_lines = (S->vals / 8) + ((S->vals % 8) ? 1 : 0);
972 long total_lines = col_lines + 2 * row_lines;
973
974 char title[32];
975 sprintf(title, "SVDLIBC v. %s", SVDVersion);
976 fprintf(file, "%-72s%-8s\n", title, "<key>");
977 fprintf(file, "%14ld%14ld%14ld%14ld%14d\n", total_lines, col_lines,
978 row_lines, row_lines, 0);
979 fprintf(file, "%-14s%14ld%14ld%14ld%14d\n", "rra", S->rows, S->cols,
980 S->vals, 0);
981 fprintf(file, "%16s%16s%16s%16s\n", "(8i)", "(8i)", "(8e)", "(8e)");
982
983 for (i = 0; i <= S->cols; i++)
984 fprintf(file, "%ld%s", S->pointr[i] + 1, (((i+1) % 8) == 0) ? "\n" : " ");
985 fprintf(file, "\n");
986 for (i = 0; i < S->vals; i++)
987 fprintf(file, "%ld%s", S->rowind[i] + 1, (((i+1) % 8) == 0) ? "\n" : " ");
988 fprintf(file, "\n");
989 for (i = 0; i < S->vals; i++)
990 fprintf(file, "%g%s", S->value[i], (((i+1) % 8) == 0) ? "\n" : " ");
991 fprintf(file, "\n");
992 }
993
994
svdLoadSparseTextFile(FILE * file)995 static SMat svdLoadSparseTextFile(FILE *file) {
996 long c, i, n, v, rows, cols, vals;
997 SMat S;
998 if (fscanf(file, " %ld %ld %ld", &rows, &cols, &vals) != 3) {
999 svd_error("svdLoadSparseTextFile: bad file format");
1000 return NULL;
1001 }
1002
1003 S = svdNewSMat(rows, cols, vals);
1004 if (!S) return NULL;
1005
1006 for (c = 0, v = 0; c < cols; c++) {
1007 if (fscanf(file, " %ld", &n) != 1) {
1008 svd_error("svdLoadSparseTextFile: bad file format");
1009 return NULL;
1010 }
1011 S->pointr[c] = v;
1012 for (i = 0; i < n; i++, v++) {
1013 if (fscanf(file, " %ld %lf", S->rowind + v, S->value + v) != 2) {
1014 svd_error("svdLoadSparseTextFile: bad file format");
1015 return NULL;
1016 }
1017 }
1018 }
1019 S->pointr[cols] = vals;
1020 return S;
1021 }
1022
svdWriteSparseTextFile(SMat S,FILE * file)1023 static void svdWriteSparseTextFile(SMat S, FILE *file) {
1024 int c, v;
1025 fprintf(file, "%ld %ld %ld\n", S->rows, S->cols, S->vals);
1026 for (c = 0, v = 0; c < S->cols; c++) {
1027 fprintf(file, "%ld\n", S->pointr[c + 1] - S->pointr[c]);
1028 for (; v < S->pointr[c+1]; v++)
1029 fprintf(file, "%ld %g\n", S->rowind[v], S->value[v]);
1030 }
1031 }
1032
1033
svdLoadSparseBinaryFile(FILE * file)1034 static SMat svdLoadSparseBinaryFile(FILE *file) {
1035 int rows, cols, vals, n, c, i, v, r, e = 0;
1036 float f;
1037 SMat S;
1038 e += svd_readBinInt(file, &rows);
1039 e += svd_readBinInt(file, &cols);
1040 e += svd_readBinInt(file, &vals);
1041 if (e) {
1042 svd_error("svdLoadSparseBinaryFile: bad file format");
1043 return NULL;
1044 }
1045
1046 S = svdNewSMat(rows, cols, vals);
1047 if (!S) return NULL;
1048
1049 for (c = 0, v = 0; c < cols; c++) {
1050 if (svd_readBinInt(file, &n)) {
1051 svd_error("svdLoadSparseBinaryFile: bad file format");
1052 return NULL;
1053 }
1054 S->pointr[c] = v;
1055 for (i = 0; i < n; i++, v++) {
1056 e += svd_readBinInt(file, &r);
1057 e += svd_readBinFloat(file, &f);
1058 if (e) {
1059 svd_error("svdLoadSparseBinaryFile: bad file format");
1060 return NULL;
1061 }
1062 S->rowind[v] = r;
1063 S->value[v] = f;
1064 }
1065 }
1066 S->pointr[cols] = vals;
1067 return S;
1068 }
1069
svdWriteSparseBinaryFile(SMat S,FILE * file)1070 static void svdWriteSparseBinaryFile(SMat S, FILE *file) {
1071 int c, v;
1072 svd_writeBinInt(file, (int) S->rows);
1073 svd_writeBinInt(file, (int) S->cols);
1074 svd_writeBinInt(file, (int) S->vals);
1075 for (c = 0, v = 0; c < S->cols; c++) {
1076 svd_writeBinInt(file, (int) (S->pointr[c + 1] - S->pointr[c]));
1077 for (; v < S->pointr[c+1]; v++) {
1078 svd_writeBinInt(file, (int) S->rowind[v]);
1079 svd_writeBinFloat(file, (float) S->value[v]);
1080 }
1081 }
1082 }
1083
1084
svdLoadDenseTextFile(FILE * file)1085 static DMat svdLoadDenseTextFile(FILE *file) {
1086 long rows, cols, i, j;
1087 DMat D;
1088 if (fscanf(file, " %ld %ld", &rows, &cols) != 2) {
1089 svd_error("svdLoadDenseTextFile: bad file format");
1090 return NULL;
1091 }
1092
1093 D = svdNewDMat(rows, cols);
1094 if (!D) return NULL;
1095
1096 for (i = 0; i < rows; i++)
1097 for (j = 0; j < cols; j++) {
1098 if (fscanf(file, " %lf", &(D->value[i][j])) != 1) {
1099 svd_error("svdLoadDenseTextFile: bad file format");
1100 return NULL;
1101 }
1102 }
1103 return D;
1104 }
1105
svdWriteDenseTextFile(DMat D,FILE * file)1106 static void svdWriteDenseTextFile(DMat D, FILE *file) {
1107 int i, j;
1108 fprintf(file, "%ld %ld\n", D->rows, D->cols);
1109 for (i = 0; i < D->rows; i++)
1110 for (j = 0; j < D->cols; j++)
1111 fprintf(file, "%g%c", D->value[i][j], (j == D->cols - 1) ? '\n' : ' ');
1112 }
1113
1114
svdLoadDenseBinaryFile(FILE * file)1115 static DMat svdLoadDenseBinaryFile(FILE *file) {
1116 int rows, cols, i, j, e = 0;
1117 float f;
1118 DMat D;
1119 e += svd_readBinInt(file, &rows);
1120 e += svd_readBinInt(file, &cols);
1121 if (e) {
1122 svd_error("svdLoadDenseBinaryFile: bad file format");
1123 return NULL;
1124 }
1125
1126 D = svdNewDMat(rows, cols);
1127 if (!D) return NULL;
1128
1129 for (i = 0; i < rows; i++)
1130 for (j = 0; j < cols; j++) {
1131 if (svd_readBinFloat(file, &f)) {
1132 svd_error("svdLoadDenseBinaryFile: bad file format");
1133 return NULL;
1134 }
1135 D->value[i][j] = f;
1136 }
1137 return D;
1138 }
1139
svdWriteDenseBinaryFile(DMat D,FILE * file)1140 static void svdWriteDenseBinaryFile(DMat D, FILE *file) {
1141 int i, j;
1142 svd_writeBinInt(file, (int) D->rows);
1143 svd_writeBinInt(file, (int) D->cols);
1144 for (i = 0; i < D->rows; i++)
1145 for (j = 0; j < D->cols; j++)
1146 svd_writeBinFloat(file, (float) D->value[i][j]);
1147 }
1148
1149
svdLoadSparseMatrix(char * filename,int format)1150 SMat svdLoadSparseMatrix(char *filename, int format) {
1151 SMat S = NULL;
1152 DMat D = NULL;
1153 FILE *file = svd_fatalReadFile(filename);
1154 switch (format) {
1155 case SVD_F_STH:
1156 S = svdLoadSparseTextHBFile(file);
1157 break;
1158 case SVD_F_ST:
1159 S = svdLoadSparseTextFile(file);
1160 break;
1161 case SVD_F_SB:
1162 S = svdLoadSparseBinaryFile(file);
1163 break;
1164 case SVD_F_DT:
1165 D = svdLoadDenseTextFile(file);
1166 break;
1167 case SVD_F_DB:
1168 D = svdLoadDenseBinaryFile(file);
1169 break;
1170 default: svd_error("svdLoadSparseMatrix: unknown format %d", format);
1171 }
1172 svd_closeFile(file);
1173 if (D) {
1174 S = svdConvertDtoS(D);
1175 svdFreeDMat(D);
1176 }
1177 return S;
1178 }
1179
svdLoadDenseMatrix(char * filename,int format)1180 DMat svdLoadDenseMatrix(char *filename, int format) {
1181 SMat S = NULL;
1182 DMat D = NULL;
1183 FILE *file = svd_fatalReadFile(filename);
1184 switch (format) {
1185 case SVD_F_STH:
1186 S = svdLoadSparseTextHBFile(file);
1187 break;
1188 case SVD_F_ST:
1189 S = svdLoadSparseTextFile(file);
1190 break;
1191 case SVD_F_SB:
1192 S = svdLoadSparseBinaryFile(file);
1193 break;
1194 case SVD_F_DT:
1195 D = svdLoadDenseTextFile(file);
1196 break;
1197 case SVD_F_DB:
1198 D = svdLoadDenseBinaryFile(file);
1199 break;
1200 default: svd_error("svdLoadSparseMatrix: unknown format %d", format);
1201 }
1202 svd_closeFile(file);
1203 if (S) {
1204 D = svdConvertStoD(S);
1205 svdFreeSMat(S);
1206 }
1207 return D;
1208 }
1209
svdWriteSparseMatrix(SMat S,char * filename,int format)1210 void svdWriteSparseMatrix(SMat S, char *filename, int format) {
1211 DMat D = NULL;
1212 FILE *file = svd_writeFile(filename, FALSE);
1213 if (!file) {
1214 svd_error("svdWriteSparseMatrix: failed to write file %s\n", filename);
1215 return;
1216 }
1217 switch (format) {
1218 case SVD_F_STH:
1219 svdWriteSparseTextHBFile(S, file);
1220 break;
1221 case SVD_F_ST:
1222 svdWriteSparseTextFile(S, file);
1223 break;
1224 case SVD_F_SB:
1225 svdWriteSparseBinaryFile(S, file);
1226 break;
1227 case SVD_F_DT:
1228 D = svdConvertStoD(S);
1229 svdWriteDenseTextFile(D, file);
1230 break;
1231 case SVD_F_DB:
1232 D = svdConvertStoD(S);
1233 svdWriteDenseBinaryFile(D, file);
1234 break;
1235 default: svd_error("svdLoadSparseMatrix: unknown format %d", format);
1236 }
1237 svd_closeFile(file);
1238 if (D) svdFreeDMat(D);
1239 }
1240
svdWriteDenseMatrix(DMat D,char * filename,int format)1241 void svdWriteDenseMatrix(DMat D, char *filename, int format) {
1242 SMat S = NULL;
1243 FILE *file = svd_writeFile(filename, FALSE);
1244 if (!file) {
1245 svd_error("svdWriteDenseMatrix: failed to write file %s\n", filename);
1246 return;
1247 }
1248 switch (format) {
1249 case SVD_F_STH:
1250 S = svdConvertDtoS(D);
1251 svdWriteSparseTextHBFile(S, file);
1252 break;
1253 case SVD_F_ST:
1254 S = svdConvertDtoS(D);
1255 svdWriteSparseTextFile(S, file);
1256 break;
1257 case SVD_F_SB:
1258 S = svdConvertDtoS(D);
1259 svdWriteSparseBinaryFile(S, file);
1260 break;
1261 case SVD_F_DT:
1262 svdWriteDenseTextFile(D, file);
1263 break;
1264 case SVD_F_DB:
1265 svdWriteDenseBinaryFile(D, file);
1266 break;
1267 default: svd_error("svdLoadSparseMatrix: unknown format %d", format);
1268 }
1269 svd_closeFile(file);
1270 if (S) svdFreeSMat(S);
1271 }
1272 /*
1273 Copyright © 2002, University of Tennessee Research Foundation.
1274 All rights reserved.
1275
1276 Redistribution and use in source and binary forms, with or without
1277 modification, are permitted provided that the following conditions are met:
1278
1279 * Redistributions of source code must retain the above copyright notice, this
1280 list of conditions and the following disclaimer.
1281
1282 Redistributions in binary form must reproduce the above copyright notice,
1283 this list of conditions and the following disclaimer in the documentation
1284 and/or other materials provided with the distribution.
1285
1286 * Neither the name of the University of Tennessee nor the names of its
1287 contributors may be used to endorse or promote products derived from this
1288 software without specific prior written permission.
1289
1290 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
1291 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
1292 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
1293 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
1294 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
1295 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
1296 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
1297 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
1298 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
1299 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
1300 POSSIBILITY OF SUCH DAMAGE.
1301 */
1302
1303
1304 #define MAXLL 2
1305
1306 #define LMTNW 100000000 /* max. size of working area allowed */
1307
1308 enum storeVals {STORQ = 1, RETRQ, STORP, RETRP};
1309
1310 static char *error_msg[] = { /* error messages used by function *
1311 * check_parameters */
1312 NULL,
1313 "",
1314 "ENDL MUST BE LESS THAN ENDR",
1315 "REQUESTED DIMENSIONS CANNOT EXCEED NUM ITERATIONS",
1316 "ONE OF YOUR DIMENSIONS IS LESS THAN OR EQUAL TO ZERO",
1317 "NUM ITERATIONS (NUMBER OF LANCZOS STEPS) IS INVALID",
1318 "REQUESTED DIMENSIONS (NUMBER OF EIGENPAIRS DESIRED) IS INVALID",
1319 "6*N+4*ITERATIONS+1 + ITERATIONS*ITERATIONS CANNOT EXCEED NW",
1320 "6*N+4*ITERATIONS+1 CANNOT EXCEED NW", NULL};
1321
1322 double **LanStore, *OPBTemp;
1323 double eps, eps1, reps, eps34;
1324 long ierr;
1325 /*
1326 double rnm, anorm, tol;
1327 FILE *fp_out1, *fp_out2;
1328 */
1329
1330 void purge(long n, long ll, double *r, double *q, double *ra,
1331 double *qa, double *wrk, double *eta, double *oldeta, long step,
1332 double *rnmp, double tol);
1333 void ortbnd(double *alf, double *eta, double *oldeta, double *bet, long step,
1334 double rnm);
1335 double startv(SMat A, double *wptr[], long step, long n);
1336 void store(long, long, long, double *);
1337 void imtql2(long, long, double *, double *, double *);
1338 void imtqlb(long n, double d[], double e[], double bnd[]);
1339 void write_header(long, long, double, double, long, double, long, long,
1340 long);
1341 long check_parameters(SMat A, long dimensions, long iterations,
1342 double endl, double endr, long vectors);
1343 int lanso(SMat A, long iterations, long dimensions, double endl,
1344 double endr, double *ritz, double *bnd, double *wptr[],
1345 long *neigp, long n);
1346 long ritvec(long n, SMat A, SVDRec R, double kappa, double *ritz,
1347 double *bnd, double *alf, double *bet, double *w2,
1348 long steps, long neig);
1349 long lanczos_step(SMat A, long first, long last, double *wptr[],
1350 double *alf, double *eta, double *oldeta,
1351 double *bet, long *ll, long *enough, double *rnmp,
1352 double *tolp, long n);
1353 void stpone(SMat A, double *wrkptr[], double *rnmp, double *tolp, long n);
1354 long error_bound(long *, double, double, double *, double *, long step,
1355 double tol);
1356 void machar(long *ibeta, long *it, long *irnd, long *machep, long *negep);
1357
1358 /***********************************************************************
1359 * *
1360 * main() *
1361 * Sparse SVD(A) via Eigensystem of A'A symmetric Matrix *
1362 * (double precision) *
1363 * *
1364 ***********************************************************************/
1365 /***********************************************************************
1366
1367 Description
1368 -----------
1369
1370 This sample program uses landr to compute singular triplets of A via
1371 the equivalent symmetric eigenvalue problem
1372
1373 B x = lambda x, where x' = (u',v'), lambda = sigma**2,
1374 where sigma is a singular value of A,
1375
1376 B = A'A , and A is m (nrow) by n (ncol) (nrow >> ncol),
1377
1378 so that {u,sqrt(lambda),v} is a singular triplet of A.
1379 (A' = transpose of A)
1380
1381 User supplied routines: svd_opa, opb, store, timer
1382
1383 svd_opa( x,y) takes an n-vector x and returns A*x in y.
1384 svd_opb(ncol,x,y) takes an n-vector x and returns B*x in y.
1385
1386 Based on operation flag isw, store(n,isw,j,s) stores/retrieves
1387 to/from storage a vector of length n in s.
1388
1389 User should edit timer() with an appropriate call to an intrinsic
1390 timing routine that returns elapsed user time.
1391
1392
1393 External parameters
1394 -------------------
1395
1396 Defined and documented in las2.h
1397
1398
1399 Local parameters
1400 ----------------
1401
1402 (input)
1403 endl left end of interval containing unwanted eigenvalues of B
1404 endr right end of interval containing unwanted eigenvalues of B
1405 kappa relative accuracy of ritz values acceptable as eigenvalues
1406 of B
1407 vectors is not equal to 1
1408 r work array
1409 n dimension of the eigenproblem for matrix B (ncol)
1410 dimensions upper limit of desired number of singular triplets of A
1411 iterations upper limit of desired number of Lanczos steps
1412 nnzero number of nonzeros in A
1413 vectors 1 indicates both singular values and singular vectors are
1414 wanted and they can be found in output file lav2;
1415 0 indicates only singular values are wanted
1416
1417 (output)
1418 ritz array of ritz values
1419 bnd array of error bounds
1420 d array of singular values
1421 memory total memory allocated in bytes to solve the B-eigenproblem
1422
1423
1424 Functions used
1425 --------------
1426
1427 BLAS svd_daxpy, svd_dscal, svd_ddot
1428 USER svd_opa, svd_opb, timer
1429 MISC write_header, check_parameters
1430 LAS2 landr
1431
1432
1433 Precision
1434 ---------
1435
1436 All floating-point calculations are done in double precision;
1437 variables are declared as long and double.
1438
1439
1440 LAS2 development
1441 ----------------
1442
1443 LAS2 is a C translation of the Fortran-77 LAS2 from the SVDPACK
1444 library written by Michael W. Berry, University of Tennessee,
1445 Dept. of Computer Science, 107 Ayres Hall, Knoxville, TN, 37996-1301
1446
1447 31 Jan 1992: Date written
1448
1449 Theresa H. Do
1450 University of Tennessee
1451 Dept. of Computer Science
1452 107 Ayres Hall
1453 Knoxville, TN, 37996-1301
1454 internet: tdo@cs.utk.edu
1455
1456 ***********************************************************************/
1457
1458 /***********************************************************************
1459 * *
1460 * check_parameters() *
1461 * *
1462 ***********************************************************************/
1463 /***********************************************************************
1464
1465 Description
1466 -----------
1467 Function validates input parameters and returns error code (long)
1468
1469 Parameters
1470 ----------
1471 (input)
1472 dimensions upper limit of desired number of eigenpairs of B
1473 iterations upper limit of desired number of lanczos steps
1474 n dimension of the eigenproblem for matrix B
1475 endl left end of interval containing unwanted eigenvalues of B
1476 endr right end of interval containing unwanted eigenvalues of B
1477 vectors 1 indicates both eigenvalues and eigenvectors are wanted
1478 and they can be found in lav2; 0 indicates eigenvalues only
1479 nnzero number of nonzero elements in input matrix (matrix A)
1480
1481 ***********************************************************************/
1482
check_parameters(SMat A,long dimensions,long iterations,double endl,double endr,long vectors)1483 long check_parameters(SMat A, long dimensions, long iterations,
1484 double endl, double endr, long vectors) {
1485 long error_index;
1486 error_index = 0;
1487
1488 if (endl >/*=*/ endr) error_index = 2;
1489 else if (dimensions > iterations) error_index = 3;
1490 else if (A->cols <= 0 || A->rows <= 0) error_index = 4;
1491 /*else if (n > A->cols || n > A->rows) error_index = 1;*/
1492 else if (iterations <= 0 || iterations > A->cols || iterations > A->rows)
1493 error_index = 5;
1494 else if (dimensions <= 0 || dimensions > iterations) error_index = 6;
1495 if (error_index)
1496 svd_error("svdLAS2 parameter error: %s\n", error_msg[error_index]);
1497 return(error_index);
1498 }
1499
1500 /***********************************************************************
1501 * *
1502 * write_header() *
1503 * Function writes out header of output file containing ritz values *
1504 * *
1505 ***********************************************************************/
1506
write_header(long iterations,long dimensions,double endl,double endr,long vectors,double kappa,long nrow,long ncol,long vals)1507 void write_header(long iterations, long dimensions, double endl, double endr,
1508 long vectors, double kappa, long nrow, long ncol,
1509 long vals) {
1510 {
1511 printf("SOLVING THE [A^TA] EIGENPROBLEM\n");
1512 printf("NO. OF ROWS = %6ld\n", nrow);
1513 printf("NO. OF COLUMNS = %6ld\n", ncol);
1514 printf("NO. OF NON-ZERO VALUES = %6ld\n", vals);
1515 printf("MATRIX DENSITY = %6.2f%%\n",
1516 ((float) vals / nrow) * 100 / ncol);
1517 /* printf("ORDER OF MATRIX A = %5ld\n", n); */
1518 printf("MAX. NO. OF LANCZOS STEPS = %6ld\n", iterations);
1519 printf("MAX. NO. OF EIGENPAIRS = %6ld\n", dimensions);
1520 printf("LEFT END OF THE INTERVAL = %9.2E\n", endl);
1521 printf("RIGHT END OF THE INTERVAL = %9.2E\n", endr);
1522 printf("KAPPA = %9.2E\n", kappa);
1523 /* printf("WANT S-VECTORS? [T/F] = %c\n", (vectors) ? 'T' : 'F'); */
1524 printf("\n");
1525 }
1526 return;
1527 }
1528
1529
1530 /***********************************************************************
1531 * *
1532 * landr() *
1533 * Lanczos algorithm with selective orthogonalization *
1534 * Using Simon's Recurrence *
1535 * (double precision) *
1536 * *
1537 ***********************************************************************/
1538 /***********************************************************************
1539
1540 Description
1541 -----------
1542
1543 landr() is the LAS2 driver routine that, upon entry,
1544 (1) checks for the validity of input parameters of the
1545 B-eigenproblem
1546 (2) determines several machine constants
1547 (3) makes a Lanczos run
1548 (4) calculates B-eigenvectors (singular vectors of A) if requested
1549 by user
1550
1551
1552 arguments
1553 ---------
1554
1555 (input)
1556 n dimension of the eigenproblem for A'A
1557 iterations upper limit of desired number of Lanczos steps
1558 dimensions upper limit of desired number of eigenpairs
1559 nnzero number of nonzeros in matrix A
1560 endl left end of interval containing unwanted eigenvalues of B
1561 endr right end of interval containing unwanted eigenvalues of B
1562 vectors 1 indicates both eigenvalues and eigenvectors are wanted
1563 and they can be found in output file lav2;
1564 0 indicates only eigenvalues are wanted
1565 kappa relative accuracy of ritz values acceptable as eigenvalues
1566 of B (singular values of A)
1567 r work array
1568
1569 (output)
1570 j number of Lanczos steps actually taken
1571 neig number of ritz values stabilized
1572 ritz array to hold the ritz values
1573 bnd array to hold the error bounds
1574
1575
1576 External parameters
1577 -------------------
1578
1579 Defined and documented in las2.h
1580
1581
1582 local parameters
1583 -------------------
1584
1585 ibeta radix for the floating-point representation
1586 it number of base ibeta digits in the floating-point significand
1587 irnd floating-point addition rounded or chopped
1588 machep machine relative precision or round-off error
1589 negeps largest negative integer
1590 wptr array of pointers each pointing to a work space
1591
1592
1593 Functions used
1594 --------------
1595
1596 MISC svd_dmax, machar, check_parameters
1597 LAS2 ritvec, lanso
1598
1599 ***********************************************************************/
1600
svdLAS2A(SMat A,long dimensions)1601 SVDRec svdLAS2A(SMat A, long dimensions) {
1602 double end[2] = {-1.0e-30, 1.0e-30};
1603 double kappa = 1e-6;
1604 if (!A) {
1605 svd_error("svdLAS2A called with NULL array\n");
1606 return NULL;
1607 }
1608 return svdLAS2(A, dimensions, 0, end, kappa);
1609 }
1610
1611
svdLAS2(SMat A,long dimensions,long iterations,double end[2],double kappa)1612 SVDRec svdLAS2(SMat A, long dimensions, long iterations, double end[2],
1613 double kappa) {
1614 char transpose = FALSE;
1615 long ibeta, it, irnd, machep, negep, n, i, steps, nsig, neig, m;
1616 double *wptr[10], *ritz, *bnd;
1617 SVDRec R = NULL;
1618 ierr = 0; // reset the global error flag
1619
1620 svdResetCounters();
1621 svd_random2(NULL) ; /* RWCox */
1622
1623 m = svd_imin(A->rows, A->cols);
1624 if (dimensions <= 0 || dimensions > m)
1625 dimensions = m;
1626 if (iterations <= 0 || iterations > m)
1627 iterations = m;
1628 if (iterations < dimensions) iterations = dimensions;
1629
1630 /* Write output header */
1631 if (SVDVerbosity > 0)
1632 write_header(iterations, dimensions, end[0], end[1], TRUE, kappa, A->rows,
1633 A->cols, A->vals);
1634
1635 /* Check parameters */
1636 if (check_parameters(A, dimensions, iterations, end[0], end[1], TRUE))
1637 return NULL;
1638
1639 /* If A is wide, the SVD is computed on its transpose for speed. */
1640 if (A->cols >= A->rows * 1.2) {
1641 if (SVDVerbosity > 0) printf("TRANSPOSING THE MATRIX FOR SPEED\n");
1642 transpose = TRUE;
1643 A = svdTransposeS(A);
1644 }
1645
1646 n = A->cols;
1647 /* Compute machine precision */
1648 machar(&ibeta, &it, &irnd, &machep, &negep);
1649 eps1 = eps * sqrt((double) n);
1650 reps = sqrt(eps);
1651 eps34 = reps * sqrt(reps);
1652
1653 /* Allocate temporary space. */
1654 if (!(wptr[0] = svd_doubleArray(n, TRUE, "las2: wptr[0]"))) goto abort;
1655 if (!(wptr[1] = svd_doubleArray(n, FALSE, "las2: wptr[1]"))) goto abort;
1656 if (!(wptr[2] = svd_doubleArray(n, FALSE, "las2: wptr[2]"))) goto abort;
1657 if (!(wptr[3] = svd_doubleArray(n, FALSE, "las2: wptr[3]"))) goto abort;
1658 if (!(wptr[4] = svd_doubleArray(n, FALSE, "las2: wptr[4]"))) goto abort;
1659 if (!(wptr[5] = svd_doubleArray(n, FALSE, "las2: wptr[5]"))) goto abort;
1660 if (!(wptr[6] = svd_doubleArray(iterations, FALSE, "las2: wptr[6]")))
1661 goto abort;
1662 if (!(wptr[7] = svd_doubleArray(iterations, FALSE, "las2: wptr[7]")))
1663 goto abort;
1664 if (!(wptr[8] = svd_doubleArray(iterations, FALSE, "las2: wptr[8]")))
1665 goto abort;
1666 if (!(wptr[9] = svd_doubleArray(iterations + 1, FALSE, "las2: wptr[9]")))
1667 goto abort;
1668 /* Calloc may be unnecessary: */
1669 if (!(ritz = svd_doubleArray(iterations + 1, TRUE, "las2: ritz")))
1670 goto abort;
1671 /* Calloc may be unnecessary: */
1672 if (!(bnd = svd_doubleArray(iterations + 1, TRUE, "las2: bnd")))
1673 goto abort;
1674 memset(bnd, 127, (iterations + 1) * sizeof(double));
1675
1676 if (!(LanStore = (double **) calloc(iterations + MAXLL, sizeof(double *))))
1677 goto abort;
1678 if (!(OPBTemp = svd_doubleArray(A->rows, FALSE, "las2: OPBTemp")))
1679 goto abort;
1680
1681 /* Actually run the lanczos thing: */
1682 steps = lanso(A, iterations, dimensions, end[0], end[1], ritz, bnd, wptr,
1683 &neig, n);
1684
1685 /* Print some stuff. */
1686 if (SVDVerbosity > 0) {
1687 printf("NUMBER OF LANCZOS STEPS = %6ld\n"
1688 "RITZ VALUES STABILIZED = %6ld\n", steps + 1, neig);
1689 }
1690 if (SVDVerbosity > 2) {
1691 printf("COMPUTED RITZ VALUES (ERROR BNDS)\n");
1692 for (i = 0; i <= steps; i++)
1693 printf("# %3ld %22.14E (%11.2E) ", i + 1, ritz[i], bnd[i]);
1694 printf("\n") ;
1695 }
1696
1697 SAFE_FREE(wptr[0]);
1698 SAFE_FREE(wptr[1]);
1699 SAFE_FREE(wptr[2]);
1700 SAFE_FREE(wptr[3]);
1701 SAFE_FREE(wptr[4]);
1702 SAFE_FREE(wptr[7]);
1703 SAFE_FREE(wptr[8]);
1704
1705 /* Compute eigenvectors */
1706 kappa = svd_dmax(fabs(kappa), eps34);
1707
1708 R = svdNewSVDRec();
1709 if (!R) {
1710 svd_error("svdLAS2: allocation of R failed");
1711 goto cleanup;
1712 }
1713 R->d = /*svd_imin(nsig, dimensions)*/dimensions;
1714 R->Ut = svdNewDMat(R->d, A->rows);
1715 R->S = svd_doubleArray(R->d, TRUE, "las2: R->s");
1716 R->Vt = svdNewDMat(R->d, A->cols);
1717 if (!R->Ut || !R->S || !R->Vt) {
1718 svd_error("svdLAS2: allocation of R failed");
1719 goto cleanup;
1720 }
1721
1722 nsig = ritvec(n, A, R, kappa, ritz, bnd, wptr[6], wptr[9], wptr[5], steps,
1723 neig);
1724
1725 if (SVDVerbosity > 1) {
1726 printf("\nSINGULAR VALUES: ");
1727 svdWriteDenseArray(R->S, R->d, "-", FALSE);
1728
1729 if (SVDVerbosity > 2) {
1730 printf("\nLEFT SINGULAR VECTORS (transpose of U): ");
1731 svdWriteDenseMatrix(R->Ut, "-", SVD_F_DT);
1732
1733 printf("\nRIGHT SINGULAR VECTORS (transpose of V): ");
1734 svdWriteDenseMatrix(R->Vt, "-", SVD_F_DT);
1735 }
1736 }
1737 if (SVDVerbosity > 0) {
1738 printf("SINGULAR VALUES FOUND = %6d\n"
1739 "SIGNIFICANT VALUES = %6ld\n", R->d, nsig);
1740 }
1741
1742 cleanup:
1743 for (i = 0; i <= 9; i++)
1744 SAFE_FREE(wptr[i]);
1745 SAFE_FREE(ritz);
1746 SAFE_FREE(bnd);
1747 if (LanStore) {
1748 for (i = 0; i < iterations + MAXLL; i++)
1749 SAFE_FREE(LanStore[i]);
1750 SAFE_FREE(LanStore);
1751 }
1752 SAFE_FREE(OPBTemp);
1753
1754 /* This swaps and transposes the singular matrices if A was transposed. */
1755 if (R && transpose) {
1756 DMat T;
1757 svdFreeSMat(A);
1758 T = R->Ut;
1759 R->Ut = R->Vt;
1760 R->Vt = T;
1761 }
1762
1763 return R;
1764 abort:
1765 svd_error("svdLAS2: fatal error, aborting");
1766 return NULL;
1767 }
1768
1769
1770 /***********************************************************************
1771 * *
1772 * ritvec() *
1773 * Function computes the singular vectors of matrix A *
1774 * *
1775 ***********************************************************************/
1776 /***********************************************************************
1777
1778 Description
1779 -----------
1780
1781 This function is invoked by landr() only if eigenvectors of the A'A
1782 eigenproblem are desired. When called, ritvec() computes the
1783 singular vectors of A and writes the result to an unformatted file.
1784
1785
1786 Parameters
1787 ----------
1788
1789 (input)
1790 nrow number of rows of A
1791 steps number of Lanczos iterations performed
1792 fp_out2 pointer to unformatted output file
1793 n dimension of matrix A
1794 kappa relative accuracy of ritz values acceptable as
1795 eigenvalues of A'A
1796 ritz array of ritz values
1797 bnd array of error bounds
1798 alf array of diagonal elements of the tridiagonal matrix T
1799 bet array of off-diagonal elements of T
1800 w1, w2 work space
1801
1802 (output)
1803 xv1 array of eigenvectors of A'A (right singular vectors of A)
1804 ierr error code
1805 0 for normal return from imtql2()
1806 k if convergence did not occur for k-th eigenvalue in
1807 imtql2()
1808 nsig number of accepted ritz values based on kappa
1809
1810 (local)
1811 s work array which is initialized to the identity matrix
1812 of order (j + 1) upon calling imtql2(). After the call,
1813 s contains the orthonormal eigenvectors of the symmetric
1814 tridiagonal matrix T
1815
1816 Functions used
1817 --------------
1818
1819 BLAS svd_dscal, svd_dcopy, svd_daxpy
1820 USER store
1821 imtql2
1822
1823 ***********************************************************************/
1824
rotateArray(double * a,int size,int x)1825 void rotateArray(double *a, int size, int x) {
1826 int i, j, n, start;
1827 double t1, t2;
1828 if (x == 0) return;
1829 j = start = 0;
1830 t1 = a[0];
1831 for (i = 0; i < size; i++) {
1832 n = (j >= x) ? j - x : j + size - x;
1833 t2 = a[n];
1834 a[n] = t1;
1835 t1 = t2;
1836 j = n;
1837 if (j == start) {
1838 start = ++j;
1839 t1 = a[j];
1840 }
1841 }
1842 }
1843
ritvec(long n,SMat A,SVDRec R,double kappa,double * ritz,double * bnd,double * alf,double * bet,double * w2,long steps,long neig)1844 long ritvec(long n, SMat A, SVDRec R, double kappa, double *ritz, double *bnd,
1845 double *alf, double *bet, double *w2, long steps, long neig) {
1846 long js, jsq, i, k, /*size,*/ id2, tmp, nsig, x;
1847 double *s, *xv2, tmp0, tmp1, xnorm, *w1 = R->Vt->value[0];
1848
1849 js = steps + 1;
1850 jsq = js * js;
1851 /*size = sizeof(double) * n;*/
1852
1853 s = svd_doubleArray(jsq, TRUE, "ritvec: s");
1854 xv2 = svd_doubleArray(n, FALSE, "ritvec: xv2");
1855
1856 /* initialize s to an identity matrix */
1857 for (i = 0; i < jsq; i+= (js+1)) s[i] = 1.0;
1858 svd_dcopy(js, alf, 1, w1, -1);
1859 svd_dcopy(steps, &bet[1], 1, &w2[1], -1);
1860
1861 /* on return from imtql2(), w1 contains eigenvalues in ascending
1862 * order and s contains the corresponding eigenvectors */
1863 imtql2(js, js, w1, w2, s);
1864
1865 /*fwrite((char *)&n, sizeof(n), 1, fp_out2);
1866 fwrite((char *)&js, sizeof(js), 1, fp_out2);
1867 fwrite((char *)&kappa, sizeof(kappa), 1, fp_out2);*/
1868 /*id = 0;*/
1869 nsig = 0;
1870
1871 if (ierr) {
1872 R->d = 0;
1873 } else {
1874 x = 0;
1875 id2 = jsq - js;
1876 for (k = 0; k < js; k++) {
1877 tmp = id2;
1878 if (bnd[k] <= kappa * fabs(ritz[k]) && k > js-neig-1) {
1879 if (--x < 0) x = R->d - 1;
1880 w1 = R->Vt->value[x];
1881 for (i = 0; i < n; i++) w1[i] = 0.0;
1882 for (i = 0; i < js; i++) {
1883 store(n, RETRQ, i, w2);
1884 svd_daxpy(n, s[tmp], w2, 1, w1, 1);
1885 tmp -= js;
1886 }
1887 /*fwrite((char *)w1, size, 1, fp_out2);*/
1888
1889 /* store the w1 vector row-wise in array xv1;
1890 * size of xv1 is (steps+1) * (nrow+ncol) elements
1891 * and each vector, even though only ncol long,
1892 * will have (nrow+ncol) elements in xv1.
1893 * It is as if xv1 is a 2-d array (steps+1) by
1894 * (nrow+ncol) and each vector occupies a row */
1895
1896 /* j is the index in the R arrays, which are sorted by high to low
1897 singular values. */
1898
1899 /*for (i = 0; i < n; i++) R->Vt->value[x]xv1[id++] = w1[i];*/
1900 /*id += nrow;*/
1901 nsig++;
1902 }
1903 id2++;
1904 }
1905 SAFE_FREE(s);
1906
1907 /* Rotate the singular vectors and values. */
1908 /* x is now the location of the highest singular value. */
1909 rotateArray(R->Vt->value[0], R->Vt->rows * R->Vt->cols,
1910 x * R->Vt->cols);
1911 R->d = svd_imin(R->d, nsig);
1912 for (x = 0; x < R->d; x++) {
1913 /* multiply by matrix B first */
1914 svd_opb(A, R->Vt->value[x], xv2, OPBTemp);
1915 tmp0 = svd_ddot(n, R->Vt->value[x], 1, xv2, 1);
1916 svd_daxpy(n, -tmp0, R->Vt->value[x], 1, xv2, 1);
1917 tmp0 = sqrt(tmp0);
1918 xnorm = sqrt(svd_ddot(n, xv2, 1, xv2, 1));
1919
1920 /* multiply by matrix A to get (scaled) left s-vector */
1921 svd_opa(A, R->Vt->value[x], R->Ut->value[x]);
1922 tmp1 = 1.0 / tmp0;
1923 svd_dscal(A->rows, tmp1, R->Ut->value[x], 1);
1924 xnorm *= tmp1;
1925 bnd[i] = xnorm;
1926 R->S[x] = tmp0;
1927 }
1928 }
1929
1930 SAFE_FREE(s);
1931 SAFE_FREE(xv2);
1932 return nsig;
1933 }
1934
1935 /*----------------------------------------------------------------------*/
1936
1937 static int vstep=0 ;
vstep_print(void)1938 static void vstep_print(void)
1939 {
1940 static char xx[10] = "0123456789" ;
1941 fprintf(stderr , "%c" , xx[vstep%10] ) ;
1942 if( vstep%10 == 9) fprintf(stderr,".") ;
1943 vstep++ ;
1944 }
1945
1946 /***********************************************************************
1947 * *
1948 * lanso() *
1949 * *
1950 ***********************************************************************/
1951 /***********************************************************************
1952
1953 Description
1954 -----------
1955
1956 Function determines when the restart of the Lanczos algorithm should
1957 occur and when it should terminate.
1958
1959 Arguments
1960 ---------
1961
1962 (input)
1963 n dimension of the eigenproblem for matrix B
1964 iterations upper limit of desired number of lanczos steps
1965 dimensions upper limit of desired number of eigenpairs
1966 endl left end of interval containing unwanted eigenvalues
1967 endr right end of interval containing unwanted eigenvalues
1968 ritz array to hold the ritz values
1969 bnd array to hold the error bounds
1970 wptr array of pointers that point to work space:
1971 wptr[0]-wptr[5] six vectors of length n
1972 wptr[6] array to hold diagonal of the tridiagonal matrix T
1973 wptr[9] array to hold off-diagonal of T
1974 wptr[7] orthogonality estimate of Lanczos vectors at
1975 step j
1976 wptr[8] orthogonality estimate of Lanczos vectors at
1977 step j-1
1978
1979 (output)
1980 j number of Lanczos steps actually taken
1981 neig number of ritz values stabilized
1982 ritz array to hold the ritz values
1983 bnd array to hold the error bounds
1984 ierr (globally declared) error flag
1985 ierr = 8192 if stpone() fails to find a starting vector
1986 ierr = k if convergence did not occur for k-th eigenvalue
1987 in imtqlb()
1988 ierr = 0 otherwise
1989
1990
1991 Functions used
1992 --------------
1993
1994 LAS stpone, error_bound, lanczos_step
1995 MISC svd_dsort2
1996 UTILITY svd_imin, svd_imax
1997
1998 ***********************************************************************/
1999
lanso(SMat A,long iterations,long dimensions,double endl,double endr,double * ritz,double * bnd,double * wptr[],long * neigp,long n)2000 int lanso(SMat A, long iterations, long dimensions, double endl,
2001 double endr, double *ritz, double *bnd, double *wptr[],
2002 long *neigp, long n) {
2003 double *alf, *eta, *oldeta, *bet, *wrk, rnm, tol;
2004 long ll, first, last, ENOUGH, id2, id3, i, l, neig, j = 0, intro = 0;
2005
2006 alf = wptr[6];
2007 eta = wptr[7];
2008 oldeta = wptr[8];
2009 bet = wptr[9];
2010 wrk = wptr[5];
2011
2012 /* take the first step */
2013 stpone(A, wptr, &rnm, &tol, n);
2014 if (!rnm || ierr) return 0;
2015 eta[0] = eps1;
2016 oldeta[0] = eps1;
2017 ll = 0;
2018 first = 1;
2019 last = svd_imin(dimensions + svd_imax(8, dimensions), iterations);
2020 ENOUGH = FALSE;
2021 /*id1 = 0;*/
2022 if( SVDVerbosity > 1 ){ fprintf(stderr,"Lanczos:"); vstep=0; }
2023 while (/*id1 < dimensions && */!ENOUGH) {
2024 if (rnm <= tol) rnm = 0.0;
2025
2026 /* the actual lanczos loop */
2027 if( SVDVerbosity > 1 ) vstep_print() ;
2028 j = lanczos_step(A, first, last, wptr, alf, eta, oldeta, bet, &ll,
2029 &ENOUGH, &rnm, &tol, n);
2030 if( SVDVerbosity > 1 ) fprintf(stderr,".") ;
2031 if (ENOUGH) j = j - 1;
2032 else j = last - 1;
2033 first = j + 1;
2034 bet[j+1] = rnm;
2035
2036 /* analyze T */
2037 l = 0;
2038 for (id2 = 0; id2 < j; id2++) {
2039 if (l > j) break;
2040 for (i = l; i <= j; i++) if (!bet[i+1]) break;
2041 if (i > j) i = j;
2042
2043 /* now i is at the end of an unreduced submatrix */
2044 svd_dcopy(i-l+1, &alf[l], 1, &ritz[l], -1);
2045 svd_dcopy(i-l, &bet[l+1], 1, &wrk[l+1], -1);
2046
2047 imtqlb(i-l+1, &ritz[l], &wrk[l], &bnd[l]);
2048
2049 if (ierr) {
2050 svd_error("svdLAS2: imtqlb failed to converge (ierr = %ld)\n", ierr);
2051 svd_error(" l = %ld i = %ld\n", l, i);
2052 for (id3 = l; id3 <= i; id3++)
2053 svd_error(" %ld %lg %lg %lg\n",
2054 id3, ritz[id3], wrk[id3], bnd[id3]);
2055 }
2056 for (id3 = l; id3 <= i; id3++)
2057 bnd[id3] = rnm * fabs(bnd[id3]);
2058 l = i + 1;
2059 }
2060 if( SVDVerbosity > 1 ) fprintf(stderr,".") ;
2061
2062 /* sort eigenvalues into increasing order */
2063 svd_dsort2((j+1) / 2, j + 1, ritz, bnd);
2064
2065 /* for (i = 0; i < iterations; i++)
2066 printf("%f ", ritz[i]);
2067 printf("\n"); */
2068
2069 /* massage error bounds for very close ritz values */
2070 neig = error_bound(&ENOUGH, endl, endr, ritz, bnd, j, tol);
2071 *neigp = neig;
2072
2073 /* should we stop? */
2074 if (neig < dimensions) {
2075 if (!neig) {
2076 last = first + 9;
2077 intro = first;
2078 } else last = first + svd_imax(3, 1 + ((j - intro) * (dimensions-neig)) /
2079 neig);
2080 last = svd_imin(last, iterations);
2081 } else ENOUGH = TRUE;
2082 ENOUGH = ENOUGH || first >= iterations;
2083 /* id1++; */
2084 /* printf("id1=%d dimen=%d first=%d\n", id1, dimensions, first); */
2085 if( SVDVerbosity > 1 ) fprintf(stderr,".") ;
2086 }
2087 store(n, STORQ, j, wptr[1]);
2088 return j;
2089 }
2090
2091
2092 /***********************************************************************
2093 * *
2094 * lanczos_step() *
2095 * *
2096 ***********************************************************************/
2097 /***********************************************************************
2098
2099 Description
2100 -----------
2101
2102 Function embodies a single Lanczos step
2103
2104 Arguments
2105 ---------
2106
2107 (input)
2108 n dimension of the eigenproblem for matrix B
2109 first start of index through loop
2110 last end of index through loop
2111 wptr array of pointers pointing to work space
2112 alf array to hold diagonal of the tridiagonal matrix T
2113 eta orthogonality estimate of Lanczos vectors at step j
2114 oldeta orthogonality estimate of Lanczos vectors at step j-1
2115 bet array to hold off-diagonal of T
2116 ll number of intitial Lanczos vectors in local orthog.
2117 (has value of 0, 1 or 2)
2118 enough stop flag
2119
2120 Functions used
2121 --------------
2122
2123 BLAS svd_ddot, svd_dscal, svd_daxpy, svd_datx, svd_dcopy
2124 USER store
2125 LAS purge, ortbnd, startv
2126 UTILITY svd_imin, svd_imax
2127
2128 ***********************************************************************/
2129
lanczos_step(SMat A,long first,long last,double * wptr[],double * alf,double * eta,double * oldeta,double * bet,long * ll,long * enough,double * rnmp,double * tolp,long n)2130 long lanczos_step(SMat A, long first, long last, double *wptr[],
2131 double *alf, double *eta, double *oldeta,
2132 double *bet, long *ll, long *enough, double *rnmp,
2133 double *tolp, long n) {
2134 double t, *mid, rnm = *rnmp, tol = *tolp, anorm;
2135 long i, j;
2136
2137 if( SVDVerbosity > 1 ) fprintf(stderr,"[%d.%d]",(int)first,(int)last) ;
2138 for (j=first; j<last; j++) {
2139 mid = wptr[2];
2140 wptr[2] = wptr[1];
2141 wptr[1] = mid;
2142 mid = wptr[3];
2143 wptr[3] = wptr[4];
2144 wptr[4] = mid;
2145
2146 store(n, STORQ, j-1, wptr[2]);
2147 if (j-1 < MAXLL) store(n, STORP, j-1, wptr[4]);
2148 bet[j] = rnm;
2149
2150
2151
2152 if( SVDVerbosity > 1 ) fprintf(stderr,"a") ;
2153 /* restart if invariant subspace is found */
2154 if (!bet[j]) {
2155 if( SVDVerbosity > 1 ) fprintf(stderr,"b") ;
2156 rnm = startv(A, wptr, j, n);
2157 if (ierr) return j;
2158 if (!rnm) *enough = TRUE;
2159 }
2160 if (*enough) {
2161 if( SVDVerbosity > 1 ) fprintf(stderr,"c") ;
2162 /* added by Doug... */
2163 /* These lines fix a bug that occurs with low-rank matrices */
2164 mid = wptr[2];
2165 wptr[2] = wptr[1];
2166 wptr[1] = mid;
2167 /* ...added by Doug */
2168 break;
2169 }
2170
2171 /* take a lanczos step */
2172 t = 1.0 / rnm;
2173 svd_datx(n, t, wptr[0], 1, wptr[1], 1);
2174 svd_dscal(n, t, wptr[3], 1);
2175 svd_opb(A, wptr[3], wptr[0], OPBTemp);
2176 svd_daxpy(n, -rnm, wptr[2], 1, wptr[0], 1);
2177 alf[j] = svd_ddot(n, wptr[0], 1, wptr[3], 1);
2178 svd_daxpy(n, -alf[j], wptr[1], 1, wptr[0], 1);
2179
2180 /* orthogonalize against initial lanczos vectors */
2181 if (j <= MAXLL && (fabs(alf[j-1]) > 4.0 * fabs(alf[j])))
2182 *ll = j;
2183 for (i=0; i < svd_imin(*ll, j-1); i++) {
2184 store(n, RETRP, i, wptr[5]);
2185 t = svd_ddot(n, wptr[5], 1, wptr[0], 1);
2186 store(n, RETRQ, i, wptr[5]);
2187 svd_daxpy(n, -t, wptr[5], 1, wptr[0], 1);
2188 eta[i] = eps1;
2189 oldeta[i] = eps1;
2190 }
2191 if( SVDVerbosity > 1 ) fprintf(stderr,"d") ;
2192
2193 /* extended local reorthogonalization */
2194 t = svd_ddot(n, wptr[0], 1, wptr[4], 1);
2195 svd_daxpy(n, -t, wptr[2], 1, wptr[0], 1);
2196 if (bet[j] > 0.0) bet[j] = bet[j] + t;
2197 t = svd_ddot(n, wptr[0], 1, wptr[3], 1);
2198 svd_daxpy(n, -t, wptr[1], 1, wptr[0], 1);
2199 alf[j] = alf[j] + t;
2200 svd_dcopy(n, wptr[0], 1, wptr[4], 1);
2201 rnm = sqrt(svd_ddot(n, wptr[0], 1, wptr[4], 1));
2202 anorm = bet[j] + fabs(alf[j]) + rnm;
2203 tol = reps * anorm;
2204
2205 /* update the orthogonality bounds */
2206 ortbnd(alf, eta, oldeta, bet, j, rnm);
2207
2208 /* restore the orthogonality state when needed */
2209 purge(n, *ll, wptr[0], wptr[1], wptr[4], wptr[3], wptr[5], eta, oldeta,
2210 j, &rnm, tol);
2211 if (rnm <= tol) rnm = 0.0;
2212 if( SVDVerbosity > 1 ) fprintf(stderr,"e") ;
2213 }
2214 *rnmp = rnm;
2215 *tolp = tol;
2216 return j;
2217 }
2218
2219 /***********************************************************************
2220 * *
2221 * ortbnd() *
2222 * *
2223 ***********************************************************************/
2224 /***********************************************************************
2225
2226 Description
2227 -----------
2228
2229 Funtion updates the eta recurrence
2230
2231 Arguments
2232 ---------
2233
2234 (input)
2235 alf array to hold diagonal of the tridiagonal matrix T
2236 eta orthogonality estimate of Lanczos vectors at step j
2237 oldeta orthogonality estimate of Lanczos vectors at step j-1
2238 bet array to hold off-diagonal of T
2239 n dimension of the eigenproblem for matrix B
2240 j dimension of T
2241 rnm norm of the next residual vector
2242 eps1 roundoff estimate for dot product of two unit vectors
2243
2244 (output)
2245 eta orthogonality estimate of Lanczos vectors at step j+1
2246 oldeta orthogonality estimate of Lanczos vectors at step j
2247
2248
2249 Functions used
2250 --------------
2251
2252 BLAS svd_dswap
2253
2254 ***********************************************************************/
2255
ortbnd(double * alf,double * eta,double * oldeta,double * bet,long step,double rnm)2256 void ortbnd(double *alf, double *eta, double *oldeta, double *bet, long step,
2257 double rnm) {
2258 long i;
2259 if (step < 1) return;
2260 if (rnm) {
2261 if (step > 1) {
2262 oldeta[0] = (bet[1] * eta[1] + (alf[0]-alf[step]) * eta[0] -
2263 bet[step] * oldeta[0]) / rnm + eps1;
2264 }
2265 for (i=1; i<=step-2; i++)
2266 oldeta[i] = (bet[i+1] * eta[i+1] + (alf[i]-alf[step]) * eta[i] +
2267 bet[i] * eta[i-1] - bet[step] * oldeta[i])/rnm + eps1;
2268 }
2269 oldeta[step-1] = eps1;
2270 svd_dswap(step, oldeta, 1, eta, 1);
2271 eta[step] = eps1;
2272 return;
2273 }
2274
2275 /***********************************************************************
2276 * *
2277 * purge() *
2278 * *
2279 ***********************************************************************/
2280 /***********************************************************************
2281
2282 Description
2283 -----------
2284
2285 Function examines the state of orthogonality between the new Lanczos
2286 vector and the previous ones to decide whether re-orthogonalization
2287 should be performed
2288
2289
2290 Arguments
2291 ---------
2292
2293 (input)
2294 n dimension of the eigenproblem for matrix B
2295 ll number of intitial Lanczos vectors in local orthog.
2296 r residual vector to become next Lanczos vector
2297 q current Lanczos vector
2298 ra previous Lanczos vector
2299 qa previous Lanczos vector
2300 wrk temporary vector to hold the previous Lanczos vector
2301 eta state of orthogonality between r and prev. Lanczos vectors
2302 oldeta state of orthogonality between q and prev. Lanczos vectors
2303 j current Lanczos step
2304
2305 (output)
2306 r residual vector orthogonalized against previous Lanczos
2307 vectors
2308 q current Lanczos vector orthogonalized against previous ones
2309
2310
2311 Functions used
2312 --------------
2313
2314 BLAS svd_daxpy, svd_dcopy, svd_idamax, svd_ddot
2315 USER store
2316
2317 ***********************************************************************/
2318
purge(long n,long ll,double * r,double * q,double * ra,double * qa,double * wrk,double * eta,double * oldeta,long step,double * rnmp,double tol)2319 void purge(long n, long ll, double *r, double *q, double *ra,
2320 double *qa, double *wrk, double *eta, double *oldeta, long step,
2321 double *rnmp, double tol) {
2322 double t, tq, tr, reps1, rnm = *rnmp;
2323 long k, iteration, flag, i;
2324
2325 if (step < ll+2) return;
2326
2327 k = svd_idamax(step - (ll+1), &eta[ll], 1) + ll;
2328 if (fabs(eta[k]) > reps) {
2329 reps1 = eps1 / reps;
2330 iteration = 0;
2331 flag = TRUE;
2332 while (iteration < 2 && flag) {
2333 if (rnm > tol) {
2334
2335 /* bring in a lanczos vector t and orthogonalize both
2336 * r and q against it */
2337 tq = 0.0;
2338 tr = 0.0;
2339 for (i = ll; i < step; i++) {
2340 store(n, RETRQ, i, wrk);
2341 t = -svd_ddot(n, qa, 1, wrk, 1);
2342 tq += fabs(t);
2343 svd_daxpy(n, t, wrk, 1, q, 1);
2344 t = -svd_ddot(n, ra, 1, wrk, 1);
2345 tr += fabs(t);
2346 svd_daxpy(n, t, wrk, 1, r, 1);
2347 }
2348 svd_dcopy(n, q, 1, qa, 1);
2349 t = -svd_ddot(n, r, 1, qa, 1);
2350 tr += fabs(t);
2351 svd_daxpy(n, t, q, 1, r, 1);
2352 svd_dcopy(n, r, 1, ra, 1);
2353 rnm = sqrt(svd_ddot(n, ra, 1, r, 1));
2354 if (tq <= reps1 && tr <= reps1 * rnm) flag = FALSE;
2355 }
2356 iteration++;
2357 }
2358 for (i = ll; i <= step; i++) {
2359 eta[i] = eps1;
2360 oldeta[i] = eps1;
2361 }
2362 }
2363 *rnmp = rnm;
2364 return;
2365 }
2366
2367
2368 /***********************************************************************
2369 * *
2370 * stpone() *
2371 * *
2372 ***********************************************************************/
2373 /***********************************************************************
2374
2375 Description
2376 -----------
2377
2378 Function performs the first step of the Lanczos algorithm. It also
2379 does a step of extended local re-orthogonalization.
2380
2381 Arguments
2382 ---------
2383
2384 (input)
2385 n dimension of the eigenproblem for matrix B
2386
2387 (output)
2388 ierr error flag
2389 wptr array of pointers that point to work space that contains
2390 wptr[0] r[j]
2391 wptr[1] q[j]
2392 wptr[2] q[j-1]
2393 wptr[3] p
2394 wptr[4] p[j-1]
2395 wptr[6] diagonal elements of matrix T
2396
2397
2398 Functions used
2399 --------------
2400
2401 BLAS svd_daxpy, svd_datx, svd_dcopy, svd_ddot, svd_dscal
2402 USER store, opb
2403 LAS startv
2404
2405 ***********************************************************************/
2406
stpone(SMat A,double * wrkptr[],double * rnmp,double * tolp,long n)2407 void stpone(SMat A, double *wrkptr[], double *rnmp, double *tolp, long n) {
2408 double t, *alf, rnm, anorm;
2409 alf = wrkptr[6];
2410
2411 /* get initial vector; default is random */
2412 rnm = startv(A, wrkptr, 0, n);
2413 if (rnm == 0.0 || ierr != 0) return;
2414
2415 /* normalize starting vector */
2416 t = 1.0 / rnm;
2417 svd_datx(n, t, wrkptr[0], 1, wrkptr[1], 1);
2418 svd_dscal(n, t, wrkptr[3], 1);
2419
2420 /* take the first step */
2421 svd_opb(A, wrkptr[3], wrkptr[0], OPBTemp);
2422 alf[0] = svd_ddot(n, wrkptr[0], 1, wrkptr[3], 1);
2423 svd_daxpy(n, -alf[0], wrkptr[1], 1, wrkptr[0], 1);
2424 t = svd_ddot(n, wrkptr[0], 1, wrkptr[3], 1);
2425 svd_daxpy(n, -t, wrkptr[1], 1, wrkptr[0], 1);
2426 alf[0] += t;
2427 svd_dcopy(n, wrkptr[0], 1, wrkptr[4], 1);
2428 rnm = sqrt(svd_ddot(n, wrkptr[0], 1, wrkptr[4], 1));
2429 anorm = rnm + fabs(alf[0]);
2430 *rnmp = rnm;
2431 *tolp = reps * anorm;
2432
2433 return;
2434 }
2435
2436 /***********************************************************************
2437 * *
2438 * startv() *
2439 * *
2440 ***********************************************************************/
2441 /***********************************************************************
2442
2443 Description
2444 -----------
2445
2446 Function delivers a starting vector in r and returns |r|; it returns
2447 zero if the range is spanned, and ierr is non-zero if no starting
2448 vector within range of operator can be found.
2449
2450 Parameters
2451 ---------
2452
2453 (input)
2454 n dimension of the eigenproblem matrix B
2455 wptr array of pointers that point to work space
2456 j starting index for a Lanczos run
2457 eps machine epsilon (relative precision)
2458
2459 (output)
2460 wptr array of pointers that point to work space that contains
2461 r[j], q[j], q[j-1], p[j], p[j-1]
2462 ierr error flag (nonzero if no starting vector can be found)
2463
2464 Functions used
2465 --------------
2466
2467 BLAS svd_ddot, svd_dcopy, svd_daxpy
2468 USER svd_opb, store
2469 MISC random
2470
2471 ***********************************************************************/
2472
startv(SMat A,double * wptr[],long step,long n)2473 double startv(SMat A, double *wptr[], long step, long n) {
2474 double rnm2, *r, t;
2475 long irand;
2476 long id, i;
2477
2478 /* get initial vector; default is random */
2479 rnm2 = svd_ddot(n, wptr[0], 1, wptr[0], 1);
2480 irand = 918273 + step;
2481 r = wptr[0];
2482 for (id = 0; id < 3; id++) {
2483 if (id > 0 || step > 0 || rnm2 == 0)
2484 for (i = 0; i < n; i++) r[i] = svd_random2(&irand);
2485 svd_dcopy(n, wptr[0], 1, wptr[3], 1);
2486
2487 /* apply operator to put r in range (essential if m singular) */
2488 svd_opb(A, wptr[3], wptr[0], OPBTemp);
2489 svd_dcopy(n, wptr[0], 1, wptr[3], 1);
2490 rnm2 = svd_ddot(n, wptr[0], 1, wptr[3], 1);
2491 if (rnm2 > 0.0) break;
2492 }
2493
2494 /* fatal error */
2495 if (rnm2 <= 0.0) {
2496 ierr = 8192;
2497 return(-1);
2498 }
2499 if (step > 0) {
2500 for (i = 0; i < step; i++) {
2501 store(n, RETRQ, i, wptr[5]);
2502 t = -svd_ddot(n, wptr[3], 1, wptr[5], 1);
2503 svd_daxpy(n, t, wptr[5], 1, wptr[0], 1);
2504 }
2505
2506 /* make sure q[step] is orthogonal to q[step-1] */
2507 t = svd_ddot(n, wptr[4], 1, wptr[0], 1);
2508 svd_daxpy(n, -t, wptr[2], 1, wptr[0], 1);
2509 svd_dcopy(n, wptr[0], 1, wptr[3], 1);
2510 t = svd_ddot(n, wptr[3], 1, wptr[0], 1);
2511 if (t <= eps * rnm2) t = 0.0;
2512 rnm2 = t;
2513 }
2514 return(sqrt(rnm2));
2515 }
2516
2517 /***********************************************************************
2518 * *
2519 * error_bound() *
2520 * *
2521 ***********************************************************************/
2522 /***********************************************************************
2523
2524 Description
2525 -----------
2526
2527 Function massages error bounds for very close ritz values by placing
2528 a gap between them. The error bounds are then refined to reflect
2529 this.
2530
2531
2532 Arguments
2533 ---------
2534
2535 (input)
2536 endl left end of interval containing unwanted eigenvalues
2537 endr right end of interval containing unwanted eigenvalues
2538 ritz array to store the ritz values
2539 bnd array to store the error bounds
2540 enough stop flag
2541
2542
2543 Functions used
2544 --------------
2545
2546 BLAS svd_idamax
2547 UTILITY svd_dmin
2548
2549 ***********************************************************************/
2550
error_bound(long * enough,double endl,double endr,double * ritz,double * bnd,long step,double tol)2551 long error_bound(long *enough, double endl, double endr,
2552 double *ritz, double *bnd, long step, double tol) {
2553 long mid, i, neig;
2554 double gapl, gap;
2555
2556 /* massage error bounds for very close ritz values */
2557 mid = svd_idamax(step + 1, bnd, 1);
2558
2559 for (i=((step+1) + (step-1)) / 2; i >= mid + 1; i -= 1)
2560 if (fabs(ritz[i-1] - ritz[i]) < eps34 * fabs(ritz[i]))
2561 if (bnd[i] > tol && bnd[i-1] > tol) {
2562 bnd[i-1] = sqrt(bnd[i] * bnd[i] + bnd[i-1] * bnd[i-1]);
2563 bnd[i] = 0.0;
2564 }
2565
2566
2567 for (i=((step+1) - (step-1)) / 2; i <= mid - 1; i +=1 )
2568 if (fabs(ritz[i+1] - ritz[i]) < eps34 * fabs(ritz[i]))
2569 if (bnd[i] > tol && bnd[i+1] > tol) {
2570 bnd[i+1] = sqrt(bnd[i] * bnd[i] + bnd[i+1] * bnd[i+1]);
2571 bnd[i] = 0.0;
2572 }
2573
2574 /* refine the error bounds */
2575 neig = 0;
2576 gapl = ritz[step] - ritz[0];
2577 for (i = 0; i <= step; i++) {
2578 gap = gapl;
2579 if (i < step) gapl = ritz[i+1] - ritz[i];
2580 gap = svd_dmin(gap, gapl);
2581 if (gap > bnd[i]) bnd[i] = bnd[i] * (bnd[i] / gap);
2582 if (bnd[i] <= 16.0 * eps * fabs(ritz[i])) {
2583 neig++;
2584 if (!*enough) *enough = endl < ritz[i] && ritz[i] < endr;
2585 }
2586 }
2587 return neig;
2588 }
2589
2590 /***********************************************************************
2591 * *
2592 * imtqlb() *
2593 * *
2594 ***********************************************************************/
2595 /***********************************************************************
2596
2597 Description
2598 -----------
2599
2600 imtqlb() is a translation of a Fortran version of the Algol
2601 procedure IMTQL1, Num. Math. 12, 377-383(1968) by Martin and
2602 Wilkinson, as modified in Num. Math. 15, 450(1970) by Dubrulle.
2603 Handbook for Auto. Comp., vol.II-Linear Algebra, 241-248(1971).
2604 See also B. T. Smith et al, Eispack Guide, Lecture Notes in
2605 Computer Science, Springer-Verlag, (1976).
2606
2607 The function finds the eigenvalues of a symmetric tridiagonal
2608 matrix by the implicit QL method.
2609
2610
2611 Arguments
2612 ---------
2613
2614 (input)
2615 n order of the symmetric tridiagonal matrix
2616 d contains the diagonal elements of the input matrix
2617 e contains the subdiagonal elements of the input matrix in its
2618 last n-1 positions. e[0] is arbitrary
2619
2620 (output)
2621 d contains the eigenvalues in ascending order. if an error
2622 exit is made, the eigenvalues are correct and ordered for
2623 indices 0,1,...ierr, but may not be the smallest eigenvalues.
2624 e has been destroyed.
2625 ierr set to zero for normal return, j if the j-th eigenvalue has
2626 not been determined after 30 iterations.
2627
2628 Functions used
2629 --------------
2630
2631 UTILITY svd_fsign
2632 MISC svd_pythag
2633
2634 ***********************************************************************/
2635
imtqlb(long n,double d[],double e[],double bnd[])2636 void imtqlb(long n, double d[], double e[], double bnd[])
2637
2638 {
2639 long last, l, m, i, iteration;
2640
2641 /* various flags */
2642 long exchange, convergence, underflow;
2643
2644 double b, test, g, r, s, c, p, f;
2645
2646 if (n == 1) return;
2647 ierr = 0;
2648 bnd[0] = 1.0;
2649 last = n - 1;
2650 for (i = 1; i < n; i++) {
2651 bnd[i] = 0.0;
2652 e[i-1] = e[i];
2653 }
2654 e[last] = 0.0;
2655 for (l = 0; l < n; l++) {
2656 iteration = 0;
2657 while (iteration <= 30) {
2658 for (m = l; m < n; m++) {
2659 convergence = FALSE;
2660 if (m == last) break;
2661 else {
2662 test = fabs(d[m]) + fabs(d[m+1]);
2663 if (test + fabs(e[m]) == test) convergence = TRUE;
2664 }
2665 if (convergence) break;
2666 }
2667 p = d[l];
2668 f = bnd[l];
2669 if (m != l) {
2670 if (iteration == 30) {
2671 ierr = l;
2672 return;
2673 }
2674 iteration += 1;
2675 /*........ form shift ........*/
2676 g = (d[l+1] - p) / (2.0 * e[l]);
2677 r = svd_pythag(g, 1.0);
2678 g = d[m] - p + e[l] / (g + svd_fsign(r, g));
2679 s = 1.0;
2680 c = 1.0;
2681 p = 0.0;
2682 underflow = FALSE;
2683 i = m - 1;
2684 while (underflow == FALSE && i >= l) {
2685 f = s * e[i];
2686 b = c * e[i];
2687 r = svd_pythag(f, g);
2688 e[i+1] = r;
2689 if (r == 0.0) underflow = TRUE;
2690 else {
2691 s = f / r;
2692 c = g / r;
2693 g = d[i+1] - p;
2694 r = (d[i] - g) * s + 2.0 * c * b;
2695 p = s * r;
2696 d[i+1] = g + p;
2697 g = c * r - b;
2698 f = bnd[i+1];
2699 bnd[i+1] = s * bnd[i] + c * f;
2700 bnd[i] = c * bnd[i] - s * f;
2701 i--;
2702 }
2703 } /* end while (underflow != FALSE && i >= l) */
2704 /*........ recover from underflow .........*/
2705 if (underflow) {
2706 d[i+1] -= p;
2707 e[m] = 0.0;
2708 }
2709 else {
2710 d[l] -= p;
2711 e[l] = g;
2712 e[m] = 0.0;
2713 }
2714 } /* end if (m != l) */
2715 else {
2716
2717 /* order the eigenvalues */
2718 exchange = TRUE;
2719 if (l != 0) {
2720 i = l;
2721 while (i >= 1 && exchange == TRUE) {
2722 if (p < d[i-1]) {
2723 d[i] = d[i-1];
2724 bnd[i] = bnd[i-1];
2725 i--;
2726 }
2727 else exchange = FALSE;
2728 }
2729 }
2730 if (exchange) i = 0;
2731 d[i] = p;
2732 bnd[i] = f;
2733 iteration = 31;
2734 }
2735 } /* end while (iteration <= 30) */
2736 } /* end for (l=0; l<n; l++) */
2737 return;
2738 } /* end main */
2739
2740 /***********************************************************************
2741 * *
2742 * imtql2() *
2743 * *
2744 ***********************************************************************/
2745 /***********************************************************************
2746
2747 Description
2748 -----------
2749
2750 imtql2() is a translation of a Fortran version of the Algol
2751 procedure IMTQL2, Num. Math. 12, 377-383(1968) by Martin and
2752 Wilkinson, as modified in Num. Math. 15, 450(1970) by Dubrulle.
2753 Handbook for Auto. Comp., vol.II-Linear Algebra, 241-248(1971).
2754 See also B. T. Smith et al, Eispack Guide, Lecture Notes in
2755 Computer Science, Springer-Verlag, (1976).
2756
2757 This function finds the eigenvalues and eigenvectors of a symmetric
2758 tridiagonal matrix by the implicit QL method.
2759
2760
2761 Arguments
2762 ---------
2763
2764 (input)
2765 nm row dimension of the symmetric tridiagonal matrix
2766 n order of the matrix
2767 d contains the diagonal elements of the input matrix
2768 e contains the subdiagonal elements of the input matrix in its
2769 last n-1 positions. e[0] is arbitrary
2770 z contains the identity matrix
2771
2772 (output)
2773 d contains the eigenvalues in ascending order. if an error
2774 exit is made, the eigenvalues are correct but unordered for
2775 for indices 0,1,...,ierr.
2776 e has been destroyed.
2777 z contains orthonormal eigenvectors of the symmetric
2778 tridiagonal (or full) matrix. if an error exit is made,
2779 z contains the eigenvectors associated with the stored
2780 eigenvalues.
2781 ierr set to zero for normal return, j if the j-th eigenvalue has
2782 not been determined after 30 iterations.
2783
2784
2785 Functions used
2786 --------------
2787 UTILITY svd_fsign
2788 MISC svd_pythag
2789
2790 ***********************************************************************/
2791
imtql2(long nm,long n,double d[],double e[],double z[])2792 void imtql2(long nm, long n, double d[], double e[], double z[])
2793
2794 {
2795 long index, nnm, j, last, l, m, i, k, iteration, convergence, underflow;
2796 double b, test, g, r, s, c, p, f;
2797 if (n == 1) return;
2798 ierr = 0;
2799 last = n - 1;
2800 for (i = 1; i < n; i++) e[i-1] = e[i];
2801 e[last] = 0.0;
2802 nnm = n * nm;
2803 for (l = 0; l < n; l++) {
2804 iteration = 0;
2805
2806 /* look for small sub-diagonal element */
2807 while (iteration <= 30) {
2808 for (m = l; m < n; m++) {
2809 convergence = FALSE;
2810 if (m == last) break;
2811 else {
2812 test = fabs(d[m]) + fabs(d[m+1]);
2813 if (test + fabs(e[m]) == test) convergence = TRUE;
2814 }
2815 if (convergence) break;
2816 }
2817 if (m != l) {
2818
2819 /* set error -- no convergence to an eigenvalue after
2820 * 30 iterations. */
2821 if (iteration == 30) {
2822 ierr = l;
2823 return;
2824 }
2825 p = d[l];
2826 iteration += 1;
2827
2828 /* form shift */
2829 g = (d[l+1] - p) / (2.0 * e[l]);
2830 r = svd_pythag(g, 1.0);
2831 g = d[m] - p + e[l] / (g + svd_fsign(r, g));
2832 s = 1.0;
2833 c = 1.0;
2834 p = 0.0;
2835 underflow = FALSE;
2836 i = m - 1;
2837 while (underflow == FALSE && i >= l) {
2838 f = s * e[i];
2839 b = c * e[i];
2840 r = svd_pythag(f, g);
2841 e[i+1] = r;
2842 if (r == 0.0) underflow = TRUE;
2843 else {
2844 s = f / r;
2845 c = g / r;
2846 g = d[i+1] - p;
2847 r = (d[i] - g) * s + 2.0 * c * b;
2848 p = s * r;
2849 d[i+1] = g + p;
2850 g = c * r - b;
2851
2852 /* form vector */
2853 for (k = 0; k < nnm; k += n) {
2854 index = k + i;
2855 f = z[index+1];
2856 z[index+1] = s * z[index] + c * f;
2857 z[index] = c * z[index] - s * f;
2858 }
2859 i--;
2860 }
2861 } /* end while (underflow != FALSE && i >= l) */
2862 /*........ recover from underflow .........*/
2863 if (underflow) {
2864 d[i+1] -= p;
2865 e[m] = 0.0;
2866 }
2867 else {
2868 d[l] -= p;
2869 e[l] = g;
2870 e[m] = 0.0;
2871 }
2872 }
2873 else break;
2874 } /*...... end while (iteration <= 30) .........*/
2875 } /*...... end for (l=0; l<n; l++) .............*/
2876
2877 /* order the eigenvalues */
2878 for (l = 1; l < n; l++) {
2879 i = l - 1;
2880 k = i;
2881 p = d[i];
2882 for (j = l; j < n; j++) {
2883 if (d[j] < p) {
2884 k = j;
2885 p = d[j];
2886 }
2887 }
2888 /* ...and corresponding eigenvectors */
2889 if (k != i) {
2890 d[k] = d[i];
2891 d[i] = p;
2892 for (j = 0; j < nnm; j += n) {
2893 p = z[j+i];
2894 z[j+i] = z[j+k];
2895 z[j+k] = p;
2896 }
2897 }
2898 }
2899 return;
2900 } /*...... end main ............................*/
2901
2902 /***********************************************************************
2903 * *
2904 * machar() *
2905 * *
2906 ***********************************************************************/
2907 /***********************************************************************
2908
2909 Description
2910 -----------
2911
2912 This function is a partial translation of a Fortran-77 subroutine
2913 written by W. J. Cody of Argonne National Laboratory.
2914 It dynamically determines the listed machine parameters of the
2915 floating-point arithmetic. According to the documentation of
2916 the Fortran code, "the determination of the first three uses an
2917 extension of an algorithm due to M. Malcolm, ACM 15 (1972),
2918 pp. 949-951, incorporating some, but not all, of the improvements
2919 suggested by M. Gentleman and S. Marovich, CACM 17 (1974),
2920 pp. 276-277." The complete Fortran version of this translation is
2921 documented in W. J. Cody, "Machar: a Subroutine to Dynamically
2922 Determine Determine Machine Parameters," TOMS 14, December, 1988.
2923
2924
2925 Parameters reported
2926 -------------------
2927
2928 ibeta the radix for the floating-point representation
2929 it the number of base ibeta digits in the floating-point
2930 significand
2931 irnd 0 if floating-point addition chops
2932 1 if floating-point addition rounds, but not in the
2933 ieee style
2934 2 if floating-point addition rounds in the ieee style
2935 3 if floating-point addition chops, and there is
2936 partial underflow
2937 4 if floating-point addition rounds, but not in the
2938 ieee style, and there is partial underflow
2939 5 if floating-point addition rounds in the ieee style,
2940 and there is partial underflow
2941 machep the largest negative integer such that
2942 1.0+float(ibeta)**machep .ne. 1.0, except that
2943 machep is bounded below by -(it+3)
2944 negeps the largest negative integer such that
2945 1.0-float(ibeta)**negeps .ne. 1.0, except that
2946 negeps is bounded below by -(it+3)
2947
2948 ***********************************************************************/
2949
machar(long * ibeta,long * it,long * irnd,long * machep,long * negep)2950 void machar(long *ibeta, long *it, long *irnd, long *machep, long *negep) {
2951
2952 volatile double beta, betain, betah, a, b, ZERO, ONE, TWO, temp, tempa,
2953 temp1;
2954 long i, itemp;
2955
2956 ONE = (double) 1;
2957 TWO = ONE + ONE;
2958 ZERO = ONE - ONE;
2959
2960 a = ONE;
2961 temp1 = ONE;
2962 while (temp1 - ONE == ZERO) {
2963 a = a + a;
2964 temp = a + ONE;
2965 temp1 = temp - a;
2966 b += a; /* to prevent icc compiler error */
2967 }
2968 b = ONE;
2969 itemp = 0;
2970 while (itemp == 0) {
2971 b = b + b;
2972 temp = a + b;
2973 itemp = (long)(temp - a);
2974 }
2975 *ibeta = itemp;
2976 beta = (double) *ibeta;
2977
2978 *it = 0;
2979 b = ONE;
2980 temp1 = ONE;
2981 while (temp1 - ONE == ZERO) {
2982 *it = *it + 1;
2983 b = b * beta;
2984 temp = b + ONE;
2985 temp1 = temp - b;
2986 }
2987 *irnd = 0;
2988 betah = beta / TWO;
2989 temp = a + betah;
2990 if (temp - a != ZERO) *irnd = 1;
2991 tempa = a + beta;
2992 temp = tempa + betah;
2993 if ((*irnd == 0) && (temp - tempa != ZERO)) *irnd = 2;
2994
2995 *negep = *it + 3;
2996 betain = ONE / beta;
2997 a = ONE;
2998 for (i = 0; i < *negep; i++) a = a * betain;
2999 b = a;
3000 temp = ONE - a;
3001 while (temp-ONE == ZERO) {
3002 a = a * beta;
3003 *negep = *negep - 1;
3004 temp = ONE - a;
3005 }
3006 *negep = -(*negep);
3007
3008 *machep = -(*it) - 3;
3009 a = b;
3010 temp = ONE + a;
3011 while (temp - ONE == ZERO) {
3012 a = a * beta;
3013 *machep = *machep + 1;
3014 temp = ONE + a;
3015 }
3016 eps = a;
3017 return;
3018 }
3019
3020 /***********************************************************************
3021 * *
3022 * store() *
3023 * *
3024 ***********************************************************************/
3025 /***********************************************************************
3026
3027 Description
3028 -----------
3029
3030 store() is a user-supplied function which, based on the input
3031 operation flag, stores to or retrieves from memory a vector.
3032
3033
3034 Arguments
3035 ---------
3036
3037 (input)
3038 n length of vector to be stored or retrieved
3039 isw operation flag:
3040 isw = 1 request to store j-th Lanczos vector q(j)
3041 isw = 2 request to retrieve j-th Lanczos vector q(j)
3042 isw = 3 request to store q(j) for j = 0 or 1
3043 isw = 4 request to retrieve q(j) for j = 0 or 1
3044 s contains the vector to be stored for a "store" request
3045
3046 (output)
3047 s contains the vector retrieved for a "retrieve" request
3048
3049 Functions used
3050 --------------
3051
3052 BLAS svd_dcopy
3053
3054 ***********************************************************************/
3055
store(long n,long isw,long j,double * s)3056 void store(long n, long isw, long j, double *s) {
3057 /* printf("called store %ld %ld\n", isw, j); */
3058 switch(isw) {
3059 case STORQ:
3060 if (!LanStore[j + MAXLL]) {
3061 if (!(LanStore[j + MAXLL] = svd_doubleArray(n, FALSE, "LanStore[j]")))
3062 svd_fatalError("svdLAS2: failed to allocate LanStore[%d]", j + MAXLL);
3063 }
3064 svd_dcopy(n, s, 1, LanStore[j + MAXLL], 1);
3065 break;
3066 case RETRQ:
3067 if (!LanStore[j + MAXLL])
3068 svd_fatalError("svdLAS2: store (RETRQ) called on index %d (not allocated)",
3069 j + MAXLL);
3070 svd_dcopy(n, LanStore[j + MAXLL], 1, s, 1);
3071 break;
3072 case STORP:
3073 if (j >= MAXLL) {
3074 svd_error("svdLAS2: store (STORP) called with j >= MAXLL");
3075 break;
3076 }
3077 if (!LanStore[j]) {
3078 if (!(LanStore[j] = svd_doubleArray(n, FALSE, "LanStore[j]")))
3079 svd_fatalError("svdLAS2: failed to allocate LanStore[%d]", j);
3080 }
3081 svd_dcopy(n, s, 1, LanStore[j], 1);
3082 break;
3083 case RETRP:
3084 if (j >= MAXLL) {
3085 svd_error("svdLAS2: store (RETRP) called with j >= MAXLL");
3086 break;
3087 }
3088 if (!LanStore[j])
3089 svd_fatalError("svdLAS2: store (RETRP) called on index %d (not allocated)",
3090 j);
3091 svd_dcopy(n, LanStore[j], 1, s, 1);
3092 break;
3093 }
3094 return;
3095 }
3096
3097 /******************************************************************************/
3098 /*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
3099 /******************************************************************************/
3100
AFNI_svdLAS2(int m,int n,double * a,double * s,double * u,double * v)3101 void AFNI_svdLAS2( int m, int n, double *a, double *s, double *u, double *v )
3102 {
3103 DMat aD , uD , vD ; SMat aS ; SVDRec aSVD ;
3104 int ii , jj , dd ;
3105
3106 if( a == NULL || s == NULL || m < 1 || n < 1 ) return ;
3107 if( u == NULL || v == NULL ) return ;
3108
3109 #if 0
3110 fprintf(stderr,"enter AFNI_svdLAS2\n") ;
3111 #endif
3112
3113 /* copy input matrix into local dense matrix struct */
3114
3115 aD = svdNewDMat(m,n) ;
3116 for( ii=0 ; ii < m ; ii++ )
3117 for( jj=0 ; jj < n ; jj++ ) aD->value[ii][jj] = a[ii+jj*m] ;
3118
3119 /* convert dense matrix struct to sparse */
3120
3121 #if 0
3122 fprintf(stderr," convert dense to sparse\n") ;
3123 #endif
3124 aS = svdConvertDtoS(aD) ; svdFreeDMat(aD) ;
3125
3126 /* carry out the SVD */
3127
3128 #if 0
3129 fprintf(stderr," call svdLAS2A\n") ;
3130 #endif
3131 #ifdef USE_OMP
3132 if( omp_in_parallel() ){
3133 #pragma omp critical
3134 { aSVD = svdLAS2A(aS,0) ; svdFreeSMat(aS) ; }
3135 } else {
3136 aSVD = svdLAS2A(aS,0) ; svdFreeSMat(aS) ;
3137 }
3138 #else
3139 aSVD = svdLAS2A(aS,0) ; svdFreeSMat(aS) ;
3140 #endif
3141
3142 /* unpack the results into the output spaces */
3143
3144 #if 0
3145 fprintf(stderr," unpack results\n") ;
3146 #endif
3147 dd = aSVD->d ;
3148 uD = aSVD->Ut ; /* Ut has dd rows, each m long */
3149 vD = aSVD->Vt ; /* Vt has dd rows, each n long */
3150
3151 for( jj=0 ; jj < n ; jj++ ){
3152 if( jj < dd ){
3153 s[jj] = aSVD->S[jj] ;
3154 for( ii=0 ; ii < m ; ii++ ) u[ii+jj*m] = uD->value[jj][ii] ;
3155 for( ii=0 ; ii < n ; ii++ ) v[ii+jj*n] = vD->value[jj][ii] ;
3156 } else {
3157 s[jj] = 0.0 ;
3158 for( ii=0 ; ii < m ; ii++ ) u[ii+jj*m] = 0.0 ;
3159 for( ii=0 ; ii < n ; ii++ ) v[ii+jj*n] = 0.0 ;
3160 }
3161 }
3162
3163 svdFreeSVDRec(aSVD) ; return ;
3164 }
3165 #endif
3166