1c this file contains the following user-callable routines: 2c 3c 4c routine iddp_rsvd computes the SVD, to a specified precision, 5c of a matrix specified by routines for applying the matrix 6c and its transpose to arbitrary vectors. 7c This routine is randomized. 8c 9c 10ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 11c 12c 13c 14c 15 subroutine iddp_rsvd(lw,eps,m,n,matvect,p1t,p2t,p3t,p4t, 16 1 matvec,p1,p2,p3,p4,krank,iu,iv,is,w,ier) 17c 18c constructs a rank-krank SVD U Sigma V^T approximating a 19c to precision eps, where matvect is a routine which applies a^T 20c to an arbitrary vector, and matvec is a routine 21c which applies a to an arbitrary vector; U is an m x krank 22c matrix whose columns are orthonormal, V is an n x krank 23c matrix whose columns are orthonormal, and Sigma is a diagonal 24c krank x krank matrix whose entries are all nonnegative. 25c The entries of U are stored in w, starting at w(iu); 26c the entries of V are stored in w, starting at w(iv). 27c The diagonal entries of Sigma are stored in w, 28c starting at w(is). This routine uses a randomized algorithm. 29c 30c input: 31c lw -- maximum usable length (in real*8 elements) 32c of the array w 33c eps -- precision of the desired approximation 34c m -- number of rows in a 35c n -- number of columns in a 36c matvect -- routine which applies the transpose 37c of the matrix to be SVD'd 38c to an arbitrary vector; this routine must have 39c a calling sequence of the form 40c 41c matvect(m,x,n,y,p1t,p2t,p3t,p4t), 42c 43c where m is the length of x, 44c x is the vector to which the transpose 45c of the matrix is to be applied, 46c n is the length of y, 47c y is the product of the transposed matrix and x, 48c and p1t, p2t, p3t, and p4t are user-specified 49c parameters 50c p1t -- parameter to be passed to routine matvect 51c p2t -- parameter to be passed to routine matvect 52c p3t -- parameter to be passed to routine matvect 53c p4t -- parameter to be passed to routine matvect 54c matvec -- routine which applies the matrix to be SVD'd 55c to an arbitrary vector; this routine must have 56c a calling sequence of the form 57c 58c matvec(n,x,m,y,p1,p2,p3,p4), 59c 60c where n is the length of x, 61c x is the vector to which the matrix is to be applied, 62c m is the length of y, 63c y is the product of the matrix and x, 64c and p1, p2, p3, and p4 are user-specified parameters 65c p1 -- parameter to be passed to routine matvec 66c p2 -- parameter to be passed to routine matvec 67c p3 -- parameter to be passed to routine matvec 68c p4 -- parameter to be passed to routine matvec 69c 70c output: 71c krank -- rank of the SVD constructed 72c iu -- index in w of the first entry of the matrix 73c of orthonormal left singular vectors of a 74c iv -- index in w of the first entry of the matrix 75c of orthonormal right singular vectors of a 76c is -- index in w of the first entry of the array 77c of singular values of a 78c w -- array containing the singular values and singular vectors 79c of a; w doubles as a work array, and so must be at least 80c (krank+1)*(3*m+5*n+1)+25*krank**2 real*8 elements long, 81c where krank is the rank returned by the present routine 82c ier -- 0 when the routine terminates successfully; 83c -1000 when lw is too small; 84c other nonzero values when idd_id2svd bombs 85c 86c _N.B._: w must be at least (krank+1)*(3*m+5*n+1)+25*krank**2 87c real*8 elements long, where krank is the rank 88c returned by the present routine. Also, the algorithm 89c used by the present routine is randomized. 90c 91 implicit none 92 integer m,n,krank,lw,lw2,ilist,llist,iproj,icol,lcol,lp, 93 1 iwork,lwork,ier,lproj,iu,iv,is,lu,lv,ls,iui,ivi,isi,k 94 real*8 eps,p1t,p2t,p3t,p4t,p1,p2,p3,p4,w(*) 95 external matvect,matvec 96c 97c 98c Allocate some memory. 99c 100 lw2 = 0 101c 102 ilist = lw2+1 103 llist = n 104 lw2 = lw2+llist 105c 106 iproj = lw2+1 107c 108c 109c ID a. 110c 111 lp = lw-lw2 112 call iddp_rid(lp,eps,m,n,matvect,p1t,p2t,p3t,p4t,krank, 113 1 w(ilist),w(iproj),ier) 114 if(ier .ne. 0) return 115c 116c 117 if(krank .gt. 0) then 118c 119c 120c Allocate more memory. 121c 122 lproj = krank*(n-krank) 123 lw2 = lw2+lproj 124c 125 icol = lw2+1 126 lcol = m*krank 127 lw2 = lw2+lcol 128c 129 iui = lw2+1 130 lu = m*krank 131 lw2 = lw2+lu 132c 133 ivi = lw2+1 134 lv = n*krank 135 lw2 = lw2+lv 136c 137 isi = lw2+1 138 ls = krank 139 lw2 = lw2+ls 140c 141 iwork = lw2+1 142 lwork = (krank+1)*(m+3*n)+26*krank**2 143 lw2 = lw2+lwork 144c 145c 146 if(lw .lt. lw2) then 147 ier = -1000 148 return 149 endif 150c 151c 152 call iddp_rsvd0(m,n,matvect,p1t,p2t,p3t,p4t, 153 1 matvec,p1,p2,p3,p4,krank,w(iui),w(ivi), 154 2 w(isi),ier,w(ilist),w(iproj),w(icol), 155 3 w(iwork)) 156 if(ier .ne. 0) return 157c 158c 159 iu = 1 160 iv = iu+lu 161 is = iv+lv 162c 163c 164c Copy the singular values and singular vectors 165c into their proper locations. 166c 167 do k = 1,lu 168 w(iu+k-1) = w(iui+k-1) 169 enddo ! k 170c 171 do k = 1,lv 172 w(iv+k-1) = w(ivi+k-1) 173 enddo ! k 174c 175 do k = 1,ls 176 w(is+k-1) = w(isi+k-1) 177 enddo ! k 178c 179c 180 endif ! krank .gt. 0 181c 182c 183 return 184 end 185c 186c 187c 188c 189 subroutine iddp_rsvd0(m,n,matvect,p1t,p2t,p3t,p4t, 190 1 matvec,p1,p2,p3,p4,krank,u,v,s,ier, 191 2 list,proj,col,work) 192c 193c routine iddp_rsvd serves as a memory wrapper 194c for the present routine (please see routine iddp_rsvd 195c for further documentation). 196c 197 implicit none 198 integer m,n,krank,list(n),ier 199 real*8 p1t,p2t,p3t,p4t,p1,p2,p3,p4,u(m,krank),v(n,krank), 200 1 s(krank),proj(krank,n-krank),col(m*krank), 201 2 work((krank+1)*(m+3*n)+26*krank**2) 202 external matvect,matvec 203c 204c 205c Collect together the columns of a indexed by list into col. 206c 207 call idd_getcols(m,n,matvec,p1,p2,p3,p4,krank,list,col,work) 208c 209c 210c Convert the ID to an SVD. 211c 212 call idd_id2svd(m,krank,col,n,list,proj,u,v,s,ier,work) 213c 214c 215 return 216 end 217