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