ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdgesvd.f
Go to the documentation of this file.
1 
2  SUBROUTINE pdgesvd(JOBU,JOBVT,M,N,A,IA,JA,DESCA,S,U,IU,JU,DESCU,
3  + VT,IVT,JVT,DESCVT,WORK,LWORK,INFO)
4 *
5 * -- ScaLAPACK routine (version 1.7) --
6 * Univ. of Tennessee, Oak Ridge National Laboratory
7 * and Univ. of California Berkeley.
8 * Jan 2006
9 
10 *
11 * .. Scalar Arguments ..
12  CHARACTER JOBU,JOBVT
13  INTEGER IA,INFO,IU,IVT,JA,JU,JVT,LWORK,M,N
14 * ..
15 * .. Array Arguments ..
16  INTEGER DESCA(*),DESCU(*),DESCVT(*)
17  DOUBLE PRECISION A(*),U(*),VT(*),WORK(*)
18  DOUBLE PRECISION S(*)
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * PDGESVD computes the singular value decomposition (SVD) of an
25 * M-by-N matrix A, optionally computing the left and/or right
26 * singular vectors. The SVD is written as
27 *
28 * A = U * SIGMA * transpose(V)
29 *
30 * where SIGMA is an M-by-N matrix which is zero except for its
31 * min(M,N) diagonal elements, U is an M-by-M orthogonal matrix, and
32 * V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
33 * are the singular values of A and the columns of U and V are the
34 * corresponding right and left singular vectors, respectively. The
35 * singular values are returned in array S in decreasing order and
36 * only the first min(M,N) columns of U and rows of VT = V**T are
37 * computed.
38 *
39 * Notes
40 * =====
41 * Each global data object is described by an associated description
42 * vector. This vector stores the information required to establish
43 * the mapping between an object element and its corresponding process
44 * and memory location.
45 *
46 * Let A be a generic term for any 2D block cyclicly distributed array.
47 * Such a global array has an associated description vector DESCA.
48 * In the following comments, the character _ should be read as
49 * "of the global array".
50 *
51 * NOTATION STORED IN EXPLANATION
52 * --------------- -------------- --------------------------------------
53 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
54 * DTYPE_A = 1.
55 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
56 * the BLACS process grid A is distribu-
57 * ted over. The context itself is glo-
58 * bal, but the handle (the integer
59 * value) may vary.
60 * M_A (global) DESCA( M_ ) The number of rows in the global
61 * array A.
62 * N_A (global) DESCA( N_ ) The number of columns in the global
63 * array A.
64 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
65 * the rows of the array.
66 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
67 * the columns of the array.
68 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
69 * row of the array A is distributed.
70 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
71 * first column of the array A is
72 * distributed.
73 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
74 * array. LLD_A >= MAX(1,LOCr(M_A)).
75 *
76 * Let K be the number of rows or columns of a distributed matrix, and
77 * assume that its process grid has dimension r x c. LOCr( K ) denotes
78 * the number of elements of K that a process would receive if K were
79 * distributed over the r processes of its process column. Similarly,
80 * LOCc( K ) denotes the number of elements of K that a process would
81 * receive if K were distributed over the c processes of its process
82 * row. The values of LOCr() and LOCc() may be determined via a call
83 * to the ScaLAPACK tool function, NUMROC:
84 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
85 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
86 * An upper bound for these quantities may be computed by:
87 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
88 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
89 *
90 * Arguments
91 * =========
92 *
93 * MP = number of local rows in A and U
94 * NQ = number of local columns in A and VT
95 * SIZE = min( M, N )
96 * SIZEQ = number of local columns in U
97 * SIZEP = number of local rows in VT
98 *
99 * JOBU (global input) CHARACTER*1
100 * Specifies options for computing U:
101 * = 'V': the first SIZE columns of U (the left singular
102 * vectors) are returned in the array U;
103 * = 'N': no columns of U (no left singular vectors) are
104 * computed.
105 *
106 * JOBVT (global input) CHARACTER*1
107 * Specifies options for computing V**T:
108 * = 'V': the first SIZE rows of V**T (the right singular
109 * vectors) are returned in the array VT;
110 * = 'N': no rows of V**T (no right singular vectors) are
111 * computed.
112 *
113 * M (global input) INTEGER
114 * The number of rows of the input matrix A. M >= 0.
115 *
116 * N (global input) INTEGER
117 * The number of columns of the input matrix A. N >= 0.
118 *
119 * A (local input/workspace) block cyclic DOUBLE PRECISION
120 * array,
121 * global dimension (M, N), local dimension (MP, NQ)
122 * On exit, the contents of A are destroyed.
123 *
124 * IA (global input) INTEGER
125 * The row index in the global array A indicating the first
126 * row of sub( A ).
127 *
128 * JA (global input) INTEGER
129 * The column index in the global array A indicating the
130 * first column of sub( A ).
131 *
132 * DESCA (global input) INTEGER array of dimension DLEN_
133 * The array descriptor for the distributed matrix A.
134 *
135 * S (global output) DOUBLE PRECISION array, dimension SIZE
136 * The singular values of A, sorted so that S(i) >= S(i+1).
137 *
138 * U (local output) DOUBLE PRECISION array, local dimension
139 * (MP, SIZEQ), global dimension (M, SIZE)
140 * if JOBU = 'V', U contains the first min(m,n) columns of U
141 * if JOBU = 'N', U is not referenced.
142 *
143 * IU (global input) INTEGER
144 * The row index in the global array U indicating the first
145 * row of sub( U ).
146 *
147 * JU (global input) INTEGER
148 * The column index in the global array U indicating the
149 * first column of sub( U ).
150 *
151 * DESCU (global input) INTEGER array of dimension DLEN_
152 * The array descriptor for the distributed matrix U.
153 *
154 * VT (local output) DOUBLE PRECISION array, local dimension
155 * (SIZEP, NQ), global dimension (SIZE, N).
156 * If JOBVT = 'V', VT contains the first SIZE rows of
157 * V**T. If JOBVT = 'N', VT is not referenced.
158 *
159 * IVT (global input) INTEGER
160 * The row index in the global array VT indicating the first
161 * row of sub( VT ).
162 *
163 * JVT (global input) INTEGER
164 * The column index in the global array VT indicating the
165 * first column of sub( VT ).
166 *
167 * DESCVT (global input) INTEGER array of dimension DLEN_
168 * The array descriptor for the distributed matrix VT.
169 *
170 * WORK (local workspace/output) DOUBLE PRECISION array, dimension
171 * (LWORK)
172 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
173 *
174 * LWORK (local input) INTEGER
175 * The dimension of the array WORK.
176 *
177 * LWORK >= 1 + 6*SIZEB + MAX(WATOBD, WBDTOSVD),
178 *
179 * where SIZEB = MAX(M,N), and WATOBD and WBDTOSVD refer,
180 * respectively, to the workspace required to bidiagonalize
181 * the matrix A and to go from the bidiagonal matrix to the
182 * singular value decomposition U*S*VT.
183 *
184 * For WATOBD, the following holds:
185 *
186 * WATOBD = MAX(MAX(WPDLANGE,WPDGEBRD),
187 * MAX(WPDLARED2D,WP(pre)LARED1D)),
188 *
189 * where WPDLANGE, WPDLARED1D, WPDLARED2D, WPDGEBRD are the
190 * workspaces required respectively for the subprograms
191 * PDLANGE, PDLARED1D, PDLARED2D, PDGEBRD. Using the
192 * standard notation
193 *
194 * MP = NUMROC( M, MB, MYROW, DESCA( CTXT_ ), NPROW),
195 * NQ = NUMROC( N, NB, MYCOL, DESCA( LLD_ ), NPCOL),
196 *
197 * the workspaces required for the above subprograms are
198 *
199 * WPDLANGE = MP,
200 * WPDLARED1D = NQ0,
201 * WPDLARED2D = MP0,
202 * WPDGEBRD = NB*(MP + NQ + 1) + NQ,
203 *
204 * where NQ0 and MP0 refer, respectively, to the values obtained
205 * at MYCOL = 0 and MYROW = 0. In general, the upper limit for
206 * the workspace is given by a workspace required on
207 * processor (0,0):
208 *
209 * WATOBD <= NB*(MP0 + NQ0 + 1) + NQ0.
210 *
211 * In case of a homogeneous process grid this upper limit can
212 * be used as an estimate of the minimum workspace for every
213 * processor.
214 *
215 * For WBDTOSVD, the following holds:
216 *
217 * WBDTOSVD = SIZE*(WANTU*NRU + WANTVT*NCVT) +
218 * MAX(WDBDSQR,
219 * MAX(WANTU*WPDORMBRQLN, WANTVT*WPDORMBRPRT)),
220 *
221 * where
222 *
223 * 1, if left(right) singular vectors are wanted
224 * WANTU(WANTVT) =
225 * 0, otherwise
226 *
227 * and WDBDSQR, WPDORMBRQLN and WPDORMBRPRT refer respectively
228 * to the workspace required for the subprograms DBDSQR,
229 * PDORMBR(QLN), and PDORMBR(PRT), where QLN and PRT are the
230 * values of the arguments VECT, SIDE, and TRANS in the call
231 * to PDORMBR. NRU is equal to the local number of rows of
232 * the matrix U when distributed 1-dimensional "column" of
233 * processes. Analogously, NCVT is equal to the local number
234 * of columns of the matrix VT when distributed across
235 * 1-dimensional "row" of processes. Calling the LAPACK
236 * procedure DBDSQR requires
237 *
238 * WDBDSQR = MAX(1, 4*SIZE )
239 *
240 * on every processor. Finally,
241 *
242 * WPDORMBRQLN = MAX( (NB*(NB-1))/2, (SIZEQ+MP)*NB)+NB*NB,
243 * WPDORMBRPRT = MAX( (MB*(MB-1))/2, (SIZEP+NQ)*MB )+MB*MB,
244 *
245 * If LWORK = -1, then LWORK is global input and a workspace
246 * query is assumed; the routine only calculates the minimum
247 * size for the work array. The required workspace is returned
248 * as the first element of WORK and no error message is issued
249 * by PXERBLA.
250 *
251 *
252 * INFO (output) INTEGER
253 * = 0: successful exit
254 * < 0: if INFO = -i, the i-th argument had an illegal value
255 
256 * > 0: if DBDSQR did not converge
257 * If INFO = MIN(M,N) + 1, then PDGESVD has detected
258 * heterogeneity by finding that eigenvalues were not
259 * identical across the process grid. In this case, the
260 * accuracy of the results from PDGESVD cannot be
261 * guaranteed.
262 *
263 * =====================================================================
264 *
265 * The results of PDGEBRD, and therefore PDGESVD, may vary slightly
266 * from run to run with the same input data. If repeatability is an
267 * issue, call BLACS_SET with the appropriate option after defining
268 * the process grid.
269 *
270 * Alignment requirements
271 * ======================
272 *
273 * The routine PDGESVD inherits the same alignement requirement as
274 * the routine PDGEBRD, namely:
275 *
276 * The distributed submatrix sub( A ) must verify some alignment proper-
277 * ties, namely the following expressions should be true:
278 * ( MB_A.EQ.NB_A .AND. IROFFA.EQ.ICOFFA )
279 * where NB = MB_A = NB_A,
280 * IROFFA = MOD( IA-1, NB ), ICOFFA = MOD( JA-1, NB ),
281 *
282 * =====================================================================
283 *
284 *
285 * .. Parameters ..
286  INTEGER BLOCK_CYCLIC_2D,DLEN_,DTYPE_,CTXT_,M_,N_,MB_,NB_,RSRC_,
287  + csrc_,lld_,ithval
288  parameter(block_cyclic_2d=1,dlen_=9,dtype_=1,ctxt_=2,m_=3,n_=4,
289  + mb_=5,nb_=6,rsrc_=7,csrc_=8,lld_=9,ithval=10)
290  DOUBLE PRECISION ZERO,ONE
291  parameter(zero= (0.0d+0),one= (1.0d+0))
292 * ..
293 * .. Local Scalars ..
294  CHARACTER UPLO
295  INTEGER CONTEXTC,CONTEXTR,I,INDD,INDD2,INDE,INDE2,INDTAUP,INDTAUQ,
296  + indu,indv,indwork,ioffd,ioffe,iscale,j,k,ldu,ldvt,llwork,
297  + lwmin,maxim,mb,mp,mypcol,mypcolc,mypcolr,myprow,myprowc,
298  + myprowr,nb,ncvt,npcol,npcolc,npcolr,nprocs,nprow,nprowc,
299  + nprowr,nq,nru,SIZE,sizeb,sizep,sizepos,sizeq,wantu,wantvt,
300  + watobd,wbdtosvd,wdbdsqr,wpdgebrd,wpdlange,wpdormbrprt,
301  + wpdormbrqln
302  DOUBLE PRECISION ANRM,BIGNUM,EPS,RMAX,RMIN,SAFMIN,SIGMA,SMLNUM
303 * ..
304 * .. Local Arrays ..
305  INTEGER DESCTU(DLEN_),DESCTVT(DLEN_),IDUM1(3),IDUM2(3)
306  DOUBLE PRECISION C(1,1)
307 * ..
308 * .. External Functions ..
309  LOGICAL LSAME
310  INTEGER NUMROC
311  DOUBLE PRECISION PDLAMCH,PDLANGE
312  EXTERNAL lsame,numroc,pdlamch,pzlange
313 * ..
314 * .. External Subroutines ..
315  EXTERNAL blacs_get,blacs_gridexit,blacs_gridinfo,blacs_gridinit,
316  + chk1mat,dbdsqr,descinit,dgamn2d,dgamx2d,dscal,igamx2d,
317  + igebr2d,igebs2d,pchk1mat,pdgebrd,pdgemr2d,pdlared1d,
318  + pdlared2d,pdlascl,pdlaset,pdormbr,pxerbla
319 * ..
320 * .. Intrinsic Functions ..
321  INTRINSIC max,min,sqrt,dble
322 * ..
323 * .. Executable Statements ..
324 * This is just to keep ftnchek happy
325  IF (block_cyclic_2d*dtype_*lld_*mb_*m_*nb_*n_.LT.0) RETURN
326 *
327  CALL blacs_gridinfo(desca(ctxt_),nprow,npcol,myprow,mypcol)
328  iscale = 0
329  info = 0
330 *
331  IF (nprow.EQ.-1) THEN
332  info = - (800+ctxt_)
333  ELSE
334 *
335  SIZE = min(m,n)
336  sizeb = max(m,n)
337  nprocs = nprow*npcol
338  IF (m.GE.n) THEN
339  ioffd = ja - 1
340  ioffe = ia - 1
341  sizepos = 1
342  ELSE
343  ioffd = ia - 1
344  ioffe = ja - 1
345  sizepos = 3
346  END IF
347 *
348  IF (lsame(jobu,'V')) THEN
349  wantu = 1
350  ELSE
351  wantu = 0
352  END IF
353  IF (lsame(jobvt,'V')) THEN
354  wantvt = 1
355  ELSE
356  wantvt = 0
357  END IF
358 *
359  CALL chk1mat(m,3,n,4,ia,ja,desca,8,info)
360  IF (wantu.EQ.1) THEN
361  CALL chk1mat(m,3,SIZE,sizepos,iu,ju,descu,13,info)
362  END IF
363  IF (wantvt.EQ.1) THEN
364  CALL chk1mat(SIZE,sizepos,n,4,ivt,jvt,descvt,17,info)
365  END IF
366  CALL igamx2d(desca(ctxt_),'A',' ',1,1,info,1,1,1,-1,-1,0)
367 *
368  IF (info.EQ.0) THEN
369 *
370 * Set up pointers into the WORK array.
371 *
372  indd = 2
373  inde = indd + sizeb + ioffd
374  indd2 = inde + sizeb + ioffe
375  inde2 = indd2 + sizeb + ioffd
376 *
377  indtauq = inde2 + sizeb + ioffe
378  indtaup = indtauq + sizeb + ja - 1
379  indwork = indtaup + sizeb + ia - 1
380  llwork = lwork - indwork + 1
381 *
382 * Initialize contexts for "column" and "row" process matrices.
383 *
384  CALL blacs_get(desca(ctxt_),10,contextc)
385  CALL blacs_gridinit(contextc,'R',nprocs,1)
386  CALL blacs_gridinfo(contextc,nprowc,npcolc,myprowc,
387  + mypcolc)
388  CALL blacs_get(desca(ctxt_),10,contextr)
389  CALL blacs_gridinit(contextr,'R',1,nprocs)
390  CALL blacs_gridinfo(contextr,nprowr,npcolr,myprowr,
391  + mypcolr)
392 *
393 * Set local dimensions of matrices (this is for MB=NB=1).
394 *
395  nru = numroc(m,1,myprowc,0,nprocs)
396  ncvt = numroc(n,1,mypcolr,0,nprocs)
397  nb = desca(nb_)
398  mb = desca(mb_)
399  mp = numroc(m,mb,myprow,desca(rsrc_),nprow)
400  nq = numroc(n,nb,mypcol,desca(csrc_),npcol)
401  IF (wantvt.EQ.1) THEN
402  sizep = numroc(SIZE,descvt(mb_),myprow,descvt(rsrc_),
403  + nprow)
404  ELSE
405  sizep = 0
406  END IF
407  IF (wantu.EQ.1) THEN
408  sizeq = numroc(SIZE,descu(nb_),mypcol,descu(csrc_),
409  + npcol)
410  ELSE
411  sizeq = 0
412  END IF
413 *
414 * Transmit MAX(NQ0, MP0).
415 *
416  IF (myprow.EQ.0 .AND. mypcol.EQ.0) THEN
417  maxim = max(nq,mp)
418  CALL igebs2d(desca(ctxt_),'All',' ',1,1,maxim,1)
419  ELSE
420  CALL igebr2d(desca(ctxt_),'All',' ',1,1,maxim,1,0,0)
421  END IF
422 *
423  wpdlange = mp
424  wpdgebrd = nb* (mp+nq+1) + nq
425  watobd = max(max(wpdlange,wpdgebrd),maxim)
426 *
427  wdbdsqr = max(1,4*size)
428  wpdormbrqln = max((nb* (nb-1))/2, (sizeq+mp)*nb) + nb*nb
429  wpdormbrprt = max((mb* (mb-1))/2, (sizep+nq)*mb) + mb*mb
430  wbdtosvd = size* (wantu*nru+wantvt*ncvt) +
431  + max(wdbdsqr,max(wantu*wpdormbrqln,
432  + wantvt*wpdormbrprt))
433 *
434 * Finally, calculate required workspace.
435 *
436  lwmin = 1 + 6*sizeb + max(watobd,wbdtosvd)
437  work(1) = dble(lwmin)
438 *
439  IF (wantu.NE.1 .AND. .NOT. (lsame(jobu,'N'))) THEN
440  info = -1
441  ELSE IF (wantvt.NE.1 .AND. .NOT. (lsame(jobvt,'N'))) THEN
442  info = -2
443  ELSE IF (lwork.LT.lwmin .AND. lwork.NE.-1) THEN
444  info = -19
445  END IF
446 *
447  END IF
448 *
449  idum1(1) = wantu
450  idum1(2) = wantvt
451  IF (lwork.EQ.-1) THEN
452  idum1(3) = -1
453  ELSE
454  idum1(3) = 1
455  END IF
456  idum2(1) = 1
457  idum2(2) = 2
458  idum2(3) = 19
459  CALL pchk1mat(m,3,n,4,ia,ja,desca,8,3,idum1,idum2,info)
460  IF (info.EQ.0) THEN
461  IF (wantu.EQ.1) THEN
462  CALL pchk1mat(m,3,SIZE,4,iu,ju,descu,13,0,idum1,idum2,
463  + info)
464  END IF
465  IF (wantvt.EQ.1) THEN
466  CALL pchk1mat(SIZE,3,n,4,ivt,jvt,descvt,17,0,idum1,
467  + idum2,info)
468  END IF
469  END IF
470 *
471  END IF
472 *
473  IF (info.NE.0) THEN
474  CALL pxerbla(desca(ctxt_),'PDGESVD',-info)
475  RETURN
476  ELSE IF (lwork.EQ.-1) THEN
477  GO TO 40
478  END IF
479 *
480 * Quick return if possible.
481 *
482  IF (m.LE.0 .OR. n.LE.0) GO TO 40
483 *
484 * Get machine constants.
485 *
486  safmin = pdlamch(desca(ctxt_),'Safe minimum')
487  eps = pdlamch(desca(ctxt_),'Precision')
488  smlnum = safmin/eps
489  bignum = one/smlnum
490  rmin = sqrt(smlnum)
491  rmax = min(sqrt(bignum),one/sqrt(sqrt(safmin)))
492 *
493 * Scale matrix to allowable range, if necessary.
494 *
495  anrm = pdlange('1',m,n,a,ia,ja,desca,work(indwork))
496  IF (anrm.GT.zero .AND. anrm.LT.rmin) THEN
497  iscale = 1
498  sigma = rmin/anrm
499  ELSE IF (anrm.GT.rmax) THEN
500  iscale = 1
501  sigma = rmax/anrm
502  END IF
503 *
504  IF (iscale.EQ.1) THEN
505  CALL pdlascl('G',one,sigma,m,n,a,ia,ja,desca,info)
506  END IF
507 *
508  CALL pdgebrd(m,n,a,ia,ja,desca,work(indd),work(inde),
509  + work(indtauq),work(indtaup),work(indwork),llwork,
510  + info)
511 *
512 * Copy D and E to all processes.
513 * Array D is in local array of dimension:
514 * LOCc(JA+MIN(M,N)-1) if M >= N; LOCr(IA+MIN(M,N)-1) otherwise.
515 * Array E is in local array of dimension
516 * LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-2) otherwise.
517 *
518  IF (m.GE.n) THEN
519 * Distribute D
520  CALL pdlared1d(n+ioffd,ia,ja,desca,work(indd),work(indd2),
521  + work(indwork),llwork)
522 * Distribute E
523  CALL pdlared2d(m+ioffe,ia,ja,desca,work(inde),work(inde2),
524  + work(indwork),llwork)
525  ELSE
526 * Distribute D
527  CALL pdlared2d(m+ioffd,ia,ja,desca,work(indd),work(indd2),
528  + work(indwork),llwork)
529 * Distribute E
530  CALL pdlared1d(n+ioffe,ia,ja,desca,work(inde),work(inde2),
531  + work(indwork),llwork)
532  END IF
533 *
534 * Prepare for calling PDBDSQR.
535 *
536  IF (m.GE.n) THEN
537  uplo = 'U'
538  ELSE
539  uplo = 'L'
540  END IF
541 *
542  indu = indwork
543  indv = indu + size*nru*wantu
544  indwork = indv + size*ncvt*wantvt
545 *
546  ldu = max(1,nru)
547  ldvt = max(1,size)
548 *
549  CALL descinit(desctu,m,SIZE,1,1,0,0,contextc,ldu,info)
550  CALL descinit(desctvt,SIZE,n,1,1,0,0,contextr,ldvt,info)
551 *
552  IF (wantu.EQ.1) THEN
553  CALL pdlaset('Full',m,SIZE,zero,one,work(indu),1,1,desctu)
554  ELSE
555  nru = 0
556  END IF
557 *
558  IF (wantvt.EQ.1) THEN
559  CALL pdlaset('Full',SIZE,n,zero,one,work(indv),1,1,desctvt)
560  ELSE
561  ncvt = 0
562  END IF
563 *
564  CALL dbdsqr(uplo,SIZE,ncvt,nru,0,work(indd2+ioffd),
565  + work(inde2+ioffe),work(indv),SIZE,work(indu),ldu,c,1,
566  + work(indwork),info)
567 *
568 * Redistribute elements of U and VT in the block-cyclic fashion.
569 *
570  IF (wantu.EQ.1) CALL pdgemr2d(m,SIZE,work(indu),1,1,desctu,u,iu,
571  + ju,descu,descu(ctxt_))
572 *
573  IF (wantvt.EQ.1) CALL pdgemr2d(SIZE,n,work(indv),1,1,desctvt,vt,
574  + ivt,jvt,descvt,descvt(ctxt_))
575 *
576 * Set to ZERO "non-square" elements of the larger matrices U, VT.
577 *
578  IF (m.GT.n .AND. wantu.EQ.1) THEN
579  CALL pdlaset('Full',m-SIZE,SIZE,zero,zero,u,ia+SIZE,ju,descu)
580  ELSE IF (n.GT.m .AND. wantvt.EQ.1) THEN
581  CALL pdlaset('Full',SIZE,n-SIZE,zero,zero,vt,ivt,jvt+SIZE,
582  + descvt)
583  END IF
584 *
585 * Multiply Householder rotations from bidiagonalized matrix.
586 *
587  IF (wantu.EQ.1) CALL pdormbr('Q','L','N',m,SIZE,n,a,ia,ja,desca,
588  + work(indtauq),u,iu,ju,descu,
589  + work(indwork),llwork,info)
590 *
591  IF (wantvt.EQ.1) CALL pdormbr('P','R','T',SIZE,n,m,a,ia,ja,desca,
592  + work(indtaup),vt,ivt,jvt,descvt,
593  + work(indwork),llwork,info)
594 *
595 * Copy singular values into output array S.
596 *
597  DO 10 i = 1,SIZE
598  s(i) = work(indd2+ioffd+i-1)
599  10 CONTINUE
600 *
601 * If matrix was scaled, then rescale singular values appropriately.
602 *
603  IF (iscale.EQ.1) THEN
604  CALL dscal(SIZE,one/sigma,s,1)
605  END IF
606 *
607 * Compare every ith eigenvalue, or all if there are only a few,
608 * across the process grid to check for heterogeneity.
609 *
610  IF (size.LE.ithval) THEN
611  j = SIZE
612  k = 1
613  ELSE
614  j = size/ithval
615  k = ithval
616  END IF
617 *
618  DO 20 i = 1,j
619  work(i+inde) = s((i-1)*k+1)
620  work(i+indd2) = s((i-1)*k+1)
621  20 CONTINUE
622 *
623  CALL dgamn2d(desca(ctxt_),'a',' ',j,1,work(1+inde),j,1,1,-1,-1,0)
624  CALL dgamx2d(desca(ctxt_),'a',' ',j,1,work(1+indd2),j,1,1,-1,-1,0)
625 *
626  DO 30 i = 1,j
627  IF ((work(i+inde)-work(i+indd2)).NE.zero) THEN
628  info = SIZE + 1
629  END IF
630  30 CONTINUE
631 *
632  40 CONTINUE
633 *
634  CALL blacs_gridexit(contextc)
635  CALL blacs_gridexit(contextr)
636 *
637 * End of PDGESVD
638 *
639  RETURN
640  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdlared2d
subroutine pdlared2d(N, IA, JA, DESC, BYROW, BYALL, WORK, LWORK)
Definition: pdlared2d.f:2
pzlange
double precision function pzlange(NORM, M, N, A, IA, JA, DESCA, WORK)
Definition: pzlange.f:3
pdgesvd
subroutine pdgesvd(JOBU, JOBVT, M, N, A, IA, JA, DESCA, S, U, IU, JU, DESCU, VT, IVT, JVT, DESCVT, WORK, LWORK, INFO)
Definition: pdgesvd.f:4
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
pdgebrd
subroutine pdgebrd(M, N, A, IA, JA, DESCA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
Definition: pdgebrd.f:3
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
pdlascl
subroutine pdlascl(TYPE, CFROM, CTO, M, N, A, IA, JA, DESCA, INFO)
Definition: pdlascl.f:3
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pdormbr
subroutine pdormbr(VECT, SIDE, TRANS, M, N, K, A, IA, JA, DESCA, TAU, C, IC, JC, DESCC, WORK, LWORK, INFO)
Definition: pdormbr.f:3
pdlared1d
subroutine pdlared1d(N, IA, JA, DESC, BYCOL, BYALL, WORK, LWORK)
Definition: pdlared1d.f:2
min
#define min(A, B)
Definition: pcgemr.c:181