SISVD-1: SPARSE SVD VIA SUBPSPACE ITERATION USING 2-CYCLIC MATRICES SISVD IS A PROGRAM WRITTEN IN FORTRAN 77 DESIGNED TO FIND SEVERAL OF THE LARGEST EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC POSITIVE DEFINITE MATRIX B. THIS IS A MODIFIED VERSION OF THE RITZIT PROGRAM (ALGOL) ORIGINALLY DESIGNED BY RUTISHAUSER IN 1970 (NUM. MATH. 16, 205-223 AND HANDBOOK FOR AUTO. COMP., VOL.II- LINEAR ALGEBRA, 284-302), AND RECODED IN FORTRAN BY B. GARBOW (ARGONNE NATIONAL LAB) IN 1976. THE MATRIX B IS ASSUMED TO BE OF THE FORM B = [ALPHA*I A ], WHERE A IS M BY N (M>>N) AND SPARSE, [A' ALPHA*I] AND ALPHA IS AN UPPER BOUND FOR THE LARGEST SINGULAR VALUE OF THE MATRIX A. HENCE, THE SINGULAR TRIPLETS OF A ARE COMPUTED AS THE EIGENPAIRS OF B. IF SIGMA IS A SINGULAR VALUE OF THE MATRIX A, THEN (ALPHA+SIGMA) AND (ALPHA-SIGMA) ARE EIGENVALUES OF B. THE FIRST M COMPONENTS OF THE EIGENVECTORS CORRESPOND TO THE LEFT SINGULAR VECTORS OF A, AND THE LAST N COMPONENTS OF THE EIGENVECTORS OF B CORRESPOND TO THE RIGHT SINGULAR VECTORS OF A. THIS PARTICULAR IMPLEMENTATION IS DISCUSSED IN "MULTIPROCESSOR SPARSE SVD ALGORITHMS AND APPLICATIONS", PH.D. THESIS BY M. BERRY, UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN, OCTOBER 1990. - HOW TO USE SUBROUTINE RITZIT FOR THE SVD (CALLED BY MAIN PROGRAM) THE CALLING SEQUENCE FOR RITZIT IS . SUBROUTINE RITZIT(NM,N,KP,KM,EPS,OPB,INF,KEM,X,D,U,W,B,F,CX, RQ,S,IMEM) THE USER SPECIFIES AS PART OF THE PARAMETER LIST: . NM ... ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER X AS DECLARED IN MAIN PROGRAM {INTEGER}. N ... ORDER OF MATRIX B FOR SVD PROBLEM {INTEGER}. (N MUST NOT BE GREATER THAN NM) KP ... NUMBER OF SIMULTANEOUS ITERATION VECTORS {INTEGER}. KP MUST NOT BE GREATER THAN N. KM ... MAXIMUM NUMBER OF ITERATION STEPS TO BE PERFORMED {INTEGER}. IF STARTING VECTORS FOR THE ITERATION VECTORS ARE AVAILABLE, KM SHOULD BE PREFIXED WITH A MINUS SIGN. EPS ... TOLERANCE FOR ACCEPTING EIGENVECTORS {INTEGER}. OPB ... NAME OF THE SUBROUTINE THAT DEFINES THE MATRIX B. OPB IS CALLED WITH PARAMETERS (N,U,W) AND MUST COMPUTE W=BU WITHOUT ALTERING THE VECTOR U. INF ... NAME OF THE SUBROUTINE THAT MAY BE USED TO OBTAIN INFORMATION OR EXERT CONTROL DURING EXECUTION. INF IS CALLED WITH PARAMETERS (KS,KG,KH,F,M), WHERE KS IS THE NUMBER OF THE NEXT ITERATION STEP, KG IS THE NUMBER OF ALREADY ACCEPTED EIGENVECTORS, KH IS THE NUMBER OF ALREADY ACCEPTED EIGENVALUES, F IS THE ARRAY OF ERROR QUANTITIES FOR THE VECTORS OF X. AN ELEMENT OF F HAS THE VALUE 4.0 UNTIL THE CORRESPONDING EIGENVALUE HAS BEEN ACCEPTED. M IS THE DEGREE OF THE CURRENT CHEBYSHEV POLYNOMIAL. KS,KG,KH,M ARE {INTEGER}. F IS A 1-DIM. ARRAY OF LENGTH N {REAL*8}. KEM ... NUMBER OF EIGENVALUES AND CORRESPONDING EIGENVECTORS DESIRED {INTEGER}. (KEM MUST BE LESS THAN KP) X ... CONTAINS, IF KM IS NEGATIVE, THE STARTING VALUES FOR THE ITERATION VECTORS. OTHERWISE, ITS CONTENTS ARE IGNORED AND RANDOM STARTING VALUES ARE GENERATED. . RITZIT RETURNS VIA ITS PARAMETER LIST THE FOLLOWING ITEMS: . KM ... RESET TO THE MAGNITUDE OF ITS INPUT VALUE. KEM ... RESET TO THE NUMBER OF EIGENVALUES AND EIGENVECTORS ACTUALLY ACCEPTED WITHIN THE LIMIT OF KM STEPS. IMEM ... APPROXIMATE NUMBER OF BYTES NEEDED FOR THIS RUN. X ... CONTAINS IN ITS FIRST KEM COLUMNS ORTHONORMAL EIGENVECTORS OF B CORRESPONDING TO THE EIGENVALUES IN ARRAY D. THE REMAINING COLUMNS CONTAIN APPROXIMATIONS TO FURTHER EIGENVECTORS OF B (SINGULAR VECTORS OF MATRIX A). X IS A NM BY KP 2-DIM. ARRAY {REAL*8}. D ... CONTAINS IN ITS FIRST KEM POSITIONS THE ABSOLUTELY LARGEST EIGENVALUES OF B (PERTURBED SINGULAR VALUES OF A). THE REMAINING POSITIONS CONTAIN APPROXIMATIONS TO SMALLER EIGENVALUES OF B D IS A 1-DIM. ARRAY OF LENGTH KP {REAL*8}. . RITZIT REQUIRES VIA ITS PARAMETER LIST THE FOLLOWING ITEMS: . U ... TEMPORARY 1-DIM. STORAGE ARRAY OF LENGTH NM {REAL*8}. W ... TEMPORARY 2-DIM. M BY KP STORAGE ARRAY {REAL*8}. B ... TEMPORARY 2-DIM. KP BY KP STORAGE ARRAY {REAL*8}. F ... TEMPORARY 1-DIM. STORAGE ARRAY OF LENGTH KP {REAL*8}. CX ... TEMPORARY 1-DIM. STORAGE ARRAY OF LENGTH KP {REAL*8}. S ... TEMPORARY 1-DIM. STORAGE ARRAY OF LENGTH KP {REAL*8}. RQ ... TEMPORARY 1-DIM. STORAGE ARRAY OF LENGTH KP {REAL*8}. . SPARSE MATRIX-VECTOR MULTIPLICATION SUBROUTINES OPB AND OPC MUST BE SUPPLIED BY THE USER. . THE CALLING SEQUENCE FOR OPB (WAS BRIEFLY DESCRIBED ABOVE) . CALL OPB(N,X,Y) . OPB TAKES AN N-VECTOR X AND RETURNS THE N-VECTOR Y = B*X. IN OUR TEST PROGRAM, SIS1, WE EMPLOY THE HARWELL-BOEING SPARSE MATRIX FORMAT FOR ACCESSING ELEMENTS OF THE M BY N SPARSE MATRIX A AND ITS TRANSPOSE. OTHER SPARSE MATRIX FORMATS CAN BE USED, OF COURSE. . THE CALLING SEQUENCE FOR OPC IS . CALL OPC(N,X,Y) OPC TAKES AN N-VECTOR X AND RETURNS THE M-VECTOR Y=C*X, WHERE C = [B-ALPHA*I], I IS THE (M+N)-ORDER IDENTITY MATRIX, AND ALPHA IS AN UPPER BOUND FOR THE LARGEST SINGULAR OF MATRIX A, I.E., C = [ 0 A ] . [ A' 0 ] . IN ADDITION, RITZIT REQUIRES THE FOLLOWING USER-SUPPLIED PARAMETERS FOR TEMPORARY ARRAY ALLOCATIONS. . USER-SUPPLIED PARAMETERS FOR RITZIT: . NSIG ... MAXIMUM NUMBER OF SINGULAR TRIPLETS DESIRED {INTEGER}. VECTORS ... ARE SINGULAR VECTORS DESIRED ON OUTPUT? {BOOLEAN}. ALPHA ... CURRENTLY GENERATED BY MAIN PROGRAM USING MATRIX 1-NORM BUT MAY BE OVERWRITTEN BY USER. IF SET TO ZERO (AND THE MATRIX B IS SYMMETRIC INDEFINITE}, CONVERGENCE CAN BE IMPROVED AT THE COST OF CARRYING 2*KP VECTORS (VERSUS KP ITERATON VECTORS) {REAL*8}. . THE MAIN PROGRAM WILL READ THE FOLLOWING PARMS FROM THE INPUT FILE SII5: . NAME ... DATASET NAME {CHAR}. EM ... NUMBER OF DESIRED TRIPLETS {INTEGER}. NUMEXTRA ... ADDITIONAL VECTORS TO CARRY {INTEGER}. KM ... MAXIMUM NUMBER OF ITERATIONS {INTEGER}. EPS ... TOLERANCE FOR ACCEPTING SINGULAR VECTORS {INTEGER}. VECTORS ... OUTPUT SINGULAR VECTORS? {BOOLEAN}. . THIS VERSION OF SISVD IS DESIGNED TO APPROXIMATE THE KEM-LARGEST SINGULAR TRIPLETS OF A. USERS INTERESTED IN THE KEM-SMALLEST SINGULAR TRIPLETS VIA SUBSPACE ITERATION SHOULD CONSULT SISVD-2 (SIS2 CODE). . EXPLICIT CALLS TO THE UNIX TIMER "ETIME" IS MADE IN THE SIS1 SOURCE. A CALL TO ANOTHER TIMING ROUTINE (ELAPSED CPU TIME) MADE BE NEEDED. . IF THE PARAMETER "VECTORS" (READ FROM SII5) IS SET TO "TRUE", THE UNFORMATTED OUTPUT FILE SIO8 WILL CONTAIN THE APPROXIMATE SINGULAR VECTORS WRITTEN IN THE ORDER U[1], V[1], U[2], V[2], ..., U[KEM], V[KEM]. HERE U[I] AND V[I] DENOTE THE LEFT AND RIGHT SINGULAR VECTORS, RESPECTIVELY, CORRESPONDING TO THE I-TH APPROXIMATE SINGULAR VALUE. . A SAMPLE INF ROUTINE CALLED "INTROS" HAS BEEN SUPPLIED IN THE SIS1 SOURCE. THE OUTPUT FROM INTROS (CALLED BY RITZIT) IS WRITTEN TO THE FORMATTED OUTPUT FILE SIO4. - SPECIFIC REFERENCES FOR SUBSPACE ITERATION AND RITZIT RUTISHAUSER, H., SIMULTANEOUS ITERATION METHOD FOR SYMMETRIC MATRICES, NUM. MATH. L6, 205-223 (1970). (REPRINTED IN HANDBOOK FOR AUTOMATIC COMPUTATION, VOLUME II, LINEAR ALGEBRA, J. H. WILKINSON - C. REINSCH, CONTRIBUTION II/9, 284-302, SPRINGER- VERLAG, 1971. RUTISHAUSER, H., COMPUTATIONAL ASPECTS OF F. L. BAUER'S SIMULTANEOUS ITERATION METHOD., NUM. MATH. 13, 4-13 (1969). GARBOW, B. S. AND DONGARRA, J. J., PATH CHART AND DOCUMENTATION FOR THE EISPACK PACKAGE OF MATRIX EIGENSYSTEM ROUTINES, TECHNICAL MEMORANDUM NO. 250, APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY, JULY, 1974, UPDATED AUGUST, 1975. - INFORMATION PLEASE ADDRESS ALL QUESTIONS, COMMENTS, OR CORRECTIONS TO: M. W. BERRY DEPARMENT OF COMPUTER SCIENCE UNIVERSITY OF TENNESSEE 107 AYRES HALL KNOXVILLE, TN 37996-1301 EMAIL: BERRY@CS.UTK.EDU PHONE: (615) 974-5067