OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pcgesvd.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine pcgesvd (jobu, jobvt, m, n, a, ia, ja, desca, s, u, iu, ju, descu, vt, ivt, jvt, descvt, work, lwork, rwork, info)

Function/Subroutine Documentation

◆ pcgesvd()

subroutine pcgesvd ( character jobu,
character jobvt,
integer m,
integer n,
complex, dimension(*) a,
integer ia,
integer ja,
integer, dimension(*) desca,
real, dimension(*) s,
complex, dimension(*) u,
integer iu,
integer ju,
integer, dimension(*) descu,
complex, dimension(*) vt,
integer ivt,
integer jvt,
integer, dimension(*) descvt,
complex, dimension(*) work,
integer lwork,
real, dimension(*) rwork,
integer info )

Definition at line 2 of file pcgesvd.f.

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