TMS1: SPARSE SVD VIA TRACE MINIMIZATION USING 2-CYCLIC MATRICES TMS1 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. THE MATRIX B IS ASSUMED TO BE OF THE FORM B = [ALPHA*I A ], WHERE A IS NROW BY NCOL [A' ALPHA*I] (NROW>>NCOL) AND SPARSE, 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 NROW COMPONENTS OF THE EIGENVECTORS CORRESPOND TO THE LEFT SINGULAR VECTORS OF A, AND THE LAST NCOL COMPONENTS OF THE EIGENVECTORS OF B CORRESPOND TO THE RIGHT SINGULAR VECTORS OF A. A SIMILAR IMPLEMENTATION IS DISCUSSED IN "MULTIPROCESSOR SPARSE SVD ALGORITHMS AND APPLICATIONS", PH.D. THESIS BY M. BERRY, UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN, OCTOBER 1990. THIS VERSION DOES NOT IMPLEMENT CHEBYSHEV ACCELERATION. THE CODE USES RITZ-SHIFTING TO ACCELERATE CONVERGENCE. THIS IS A PARALLEL METHOD WHICH PERMITS CONCURRENT ITERATIONS OF THE CLASSICAL CONJUGATE GRADIENT METHOD. PARALLELIZATION DIRECTIVES SUCH AS THE C$DOACROSS FOR THE SEQUENT SYMMETRY MULTIPROCESSOR (SHARED MEMORY) PRECEDE LOOPS WHOSE ITERATIONS CAN BE INDEPENDENTLY SCHEDULED ACROSS PROCESSORS OF ANY PARALLEL MACHINE. THESE DIRECTIVES WILL BE TREATED AS COMMENTS BY ANY OTHER MACHINE'S PRE-PROCESSOR, OF COURSE, AND ITERATIONS OF THE CORRESPONDING LOOPS WILL BE EXECUTED IN SEQUENTIAL FASHION. - HOW TO USE SUBROUTINE TSVD1 FOR THE SVD (CALLED BY PROGRAM TMS1) THE CALLING SEQUENCE FOR TSVD1 IS . SUBROUTINE TSVD1(N,P,S,JOB,MAXI,TOL,RED,SIG,Y,LDY,MEM, TITER,ITCGT,TIME,RES,MXV,WORK1,LDW1, WORK2,LDW2,WORK3,LDW3,WORK4,LDW4, WORK5,LDW5,IWORK,LDWI,LWORK,IERR) THE USER SPECIFIES AS PART OF THE PARAMETER LIST: . N ... ORDER OF MATRIX B FOR SVD PROBLEM {INTEGER}. (N MUST NOT BE GREATER THAN SUM OF NUMBER OF ROWS AND COLUMNS OF SPARSE MATRIX A) P ... NUMBER OF DESIRED SINGULAR TRIPLETS (LARGEST) OF MATRIX A. {INTEGER}. S ... DIMENSION OF INITIAL SUBSPACE {INTEGER}. (S SHOULD BE GREATER THAN P; S=2*P IS USUALLY SAFE BUT MORE OPTIMAL RESULTS MAY BE OBTAINED IF S IS CLOSER TO P) JOB ... ACCELERATION STRATEGY SWITCH {INTEGER}. JOB = 0, NO ACCELERATION IS USED. JOB = 1, RITZ-SHIFTING IS USED. MAXI ... MAXIMUM NUMBER OF TRACE MINIMIZATION STEPS ALLOWED {INTEGER}. TOL ... USER-SUPPLIED TOLERANCE FOR RESIDUALS OF B-EIGENPAIRS WHICH APPROXIMATE A-SINGULAR TRIPLETS {DOUBLE PRECISION}. RED ... USER-SUPPLIED TOLERANCE FOR RESIDUAL REDUCTION TO INVOKE RITZ-SHIFTING (JOB=1). {DOUBLE PRECISION}. LDY ... LEADING DIMENSION OF ARRAY Y {INTEGER}. (LDY SHOULD BE GREATER OR EQUAL TO N) . TSVD1 RETURNS VIA ITS PARAMETER LIST THE FOLLOWING ITEMS: . IERR ... ERROR FLAG FOR JOB PARAMETER {INTEGER}. IERR=99, INPUT PARAMETER INVALID. IERR= 0, INPUT PARAMETER VALID. MEM ... ESTIMATE (IN BYTES) OF MEMORY REQUIRED {INTEGER}. (REAL, INTEGER, LOGICAL = 4 BYTES, DOUBLE PRECISION = 8 BYTES) MXV ... 1-DIM. ARRAY OF LENGTH 3 CONTAINING MATRIX TIMES VECTOR COUNTS {INTEGER}. MXV(1) = NUMBER OF A *X. (X IS A VECTOR) MXV(2) = NUMBER OF A'*X. MXV(3) = MXV(1) + MXV(2). SIG ... 1-DIM. ARRAY OF LENGTH P CONTAINING THE DESIRED SINGULAR VALUES OF A {DOUBLE PRECISION}. Y ... LDY BY P 2-DIM. ARRAY CONTAINING THE CORRESP. LEFT AND RIGHT SINGULAR VECTORS OF MATRIX A {DOUBLE PRECISION}. EACH COLUMN OF Y STORES THE LEFT SINGULAR VECTOR IN THE FIRST NROW ELEMENTS AND THE RIGHT SINGULAR VECTOR IN THE LAST NCOL ELEMENTS, WHERE NROW IS THE NUMBER OF ROWS OF A AND NCOL IS THE NUMBER OF COLUMNS OF A. TITER ... 1-DIM. ARRAY OF LENGTH P CONTAINING THE NUMBER OF TRACE MIN. STEPS REQUIRED FOR EACH SINGULAR TRIPLET OF A {INTEGER}. ITCGT ... 1-DIM. ARRAY OF LENGTH P CONTAINING THE NUMBER OF CONJUGATE GRADIENT STEPS TAKEN FOR EACH SINGULAR TRIPLET APPROXIMATION OF A {INTEGER}. TIME ... 1-DIM. ARRAY OF LENGTH 5 SPECIFYING TIMING BREAKDOWN (USER CPU SECONDS) {REAL}. TIME(1) = GRAM SCHMIDT ORTHOGONALIZATION. TIME(2) = SPECTRAL DECOMPOSITION. TIME(3) = CONVERGENCE CRITERIA. TIME(4) = CONJUGATE GRADIENT METHOD. TIME(5) = TOTAL TIME (SUM OF THE ABOVE). RES ... 1-DIM. ARRAY OF LENGTH P CONTAINING THE 2-NORMS OF RESIDUALS FOR THE SINGULAR TRIPLETS OF A {DOUBLE PRECISION}. . RITZIT REQUIRES VIA ITS PARAMETER LIST THE FOLLOWING ITEMS: . LDW1 ... LEADING DIMENSION OF ARRAY WORK1 {INTEGER}. (LDW1 SHOULD BE GREATER OR EQUAL TO N) WORK1 ... LDW1 BY S 2-DIM. WORK ARRAY {DOUBLE PRECISION}. LDW2 ... LEADING DIMENSION OF ARRAY WORK2 {INTEGER}. (LDW2 SHOULD BE GREATER OR EQUAL TO N) WORK2 ... LDW2 BY S 2-DIM. WORK ARRAY {DOUBLE PRECISION}. LDW3 ... LEADING DIMENSION OF ARRAY WORK3 {INTEGER}. (LDW3 SHOULD BE GREATER OR EQUAL TO S) WORK3 ... LDW3 BY S 2-DIM. WORK ARRAY {DOUBLE PRECISION}. LDW4 ... LEADING DIMENSION OF ARRAY WORK4 {INTEGER}. (LDW4 SHOULD BE GREATER OR EQUAL TO N) WORK4 ... LDW4 BY 3 2-DIM. WORK ARRAY {DOUBLE PRECISION}. LDW5 ... LEADING DIMENSION OF ARRAY WORK5 {INTEGER}. (LDW5 SHOULD BE GREATER OR EQUAL TO S) WORK5 ... LDW5 BY 4 2-DIM. WORK ARRAY {DOUBLE PRECISION}. LDWI ... LEADING DIMENSION OF ARRAY IWORK {INTEGER}. (LDWI SHOULD BE GREATER OR EQUAL TO S) IWORK ... LDWI BY 2 2-DIM. WORK ARRAY {INTEGER}. LWORK ... 1-DIM. WORK ARRAY OF LENGTH S {LOGICAL}. . SPARSE MATRIX-VECTOR MULTIPLICATION SUBROUTINES OPAT AND OPB MUST BE SUPPLIED BY THE USER. . THE CALLING SEQUENCE FOR OPAT IS . CALL OPAT(X,Y) . IF A IS AN NROW BY NCOL MATRIX, OPAT TAKES AN NROW-VECTOR, AND RETURNS THE NCOL-VECTOR Y = A'*X. IN OUR TEST PROGRAM, TMS1, WE EMPLOY THE HARWELL-BOEING SPARSE MATRIX FORMAT FOR ACCESSING ELEMENTS OF THE SPARSE MATRIX A AND ITS TRANSPOSE (A'). OTHER SPARSE MATRIX FORMATS CAN BE USED, OF COURSE. . THE CALLING SEQUENCE FOR OPB IS . CALL OPB(N,X,Y,SHIFT) OPB TAKES AN N-VECTOR X AND RETURNS THE N-VECTOR Y=D*X, WHERE D = [B-SHIFT*I], I IS THE NTH-ORDER IDENTITY MATRIX. ALPHA IS AN UPPER BOUND FOR THE LARGEST SINGULAR OF MATRIX A, HERE, N = NUMBER OF ROWS PLUS COLUMNS OF THE MATRIX A WHOSE SINGULAR TRIPLETS ARE SOUGHT. . IN ADDITION, TMS1 REQUIRES THE FOLLOWING USER-SUPPLIED PARAMETER FOR TEMPORARY ARRAY ALLOCATIONS. . USER-SUPPLIED PARAMETER FOR SUBROUTINES CGT and CGTS (BOTH CALLED BY TSVD1): . NCG ... SHOULD BE GREATER THAN OR EQUAL TO THE ORDER OF MATRIX B {INTEGER}. . THE MAIN PROGRAM (TMS1) WILL READ THE FOLLOWING PARMS FROM THE INPUT FILE TMI5: . P ... NUMBER OF SINGULAR TRIPLETS (LARGEST) {INTEGER}. S ... INITIAL SUBSPACE DIMENSION {INTEGER}. JOB ... RITZ-SHIFTING CONTROL(0=NO,1=YES) {INTEGER}. TOL ... USER-SPECIFIED RESIDUAL TOLERANCE {DOUBLE PRECISION}. RED ... RESIDUAL NORM REDUCTION FACTOR {DOUBLE PRECISION}. VEC ... OUTPUT SINGULAR VECTORS? {BOOLEAN}. . SUBROUTINE TSVD1 IS PRIMARILY DESIGNED TO APPROXIMATE THE P-LARGEST SINGULAR TRIPLETS OF A. USERS INTERESTED IN THE P-SMALLEST SINGULAR TRIPLETS VIA TRACE MINIMIZATION SHOULD USE SUBROUTINE TSVD2 (TMS2 SOURCE CODE). . EXPLICIT CALLS TO THE UNIX TIMER "DTIME" IS MADE IN THE TMS1 SOURCE. A CALL TO ANOTHER TIMING ROUTINE (DELTA CPU TIME) MADE BE NEEDED. EXPLICIT CALLS TO THE WALL-CLOCK TIMER "ELAPSE" ARE MADE IN THE MAIN PROGRAM (TMS1). THIS SUBROUTINE MAY ALSO NEED MODIFICATION. . IF THE PARAMETER "VEC" (READ FROM TMI5) IS SET TO ".TRUE.", THE UNFORMATTED OUTPUT FILE TMO8 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. . - ORIGINAL REFERENCE FOR TRACE MINIMIZATION: SAMEH, A. H. AND WISNIEWSKI, J. A., A TRACE MINIMIZATION STRATEGY FOR THE GENERALIZED EIGEVALUE PROBLEM, SIAM J. NUMER. ANAL. 19:6, 1243-1259, 1982. - 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