Routine PDGESVD() Documented By Universal Report

Routine: PDGESVD()

# lines:	639
# code:	639
# comment:	0
# blank:	0
# Variables:	90
# Callers:	0
# Callings:	9
# Words:	349
# Keywords:	210

.. Scalar Arguments ..
     ..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PDGESVD computes the singular value decomposition (SVD) of an
  M-by-N matrix A, optionally computing the left and/or right
  singular vectors. The SVD is written as
       A = U * SIGMA * transpose(V)
  where SIGMA is an M-by-N matrix which is zero except for its
  min(M,N) diagonal elements, U is an M-by-M orthogonal matrix, and
  V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA
  are the singular values of A and the columns of U and V are the
  corresponding right and left singular vectors, respectively. The
  singular values are returned in array S in decreasing order and
  only the first min(M,N) columns of U and rows of VT = V**T are
  computed.
  Notes
  =====
  Each global data object is described by an associated description
  vector. This vector stores the information required to establish
  the mapping between an object element and its corresponding process
  and memory location.
  Let A be a generic term for any 2D block cyclicly distributed array.
  Such a global array has an associated description vector DESCA.
  In the following comments, the character _ should be read as
  "of the global array".
  NOTATION        STORED IN      EXPLANATION
  --------------- -------------- --------------------------------------
  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
                                 DTYPE_A = 1.
  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
                                 the BLACS process grid A is distribu-
                                 ted over. The context itself is glo-
                                 bal, but the handle (the integer
                                 value) may vary.
  M_A    (global) DESCA( M_ )    The number of rows in the global
                                 array A.
  N_A    (global) DESCA( N_ )    The number of columns in the global
                                 array A.
  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
                                 the rows of the array.
  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
                                 the columns of the array.
  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
                                 row of the array A is distributed.
  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
                                 first column of the array A is
                                 distributed.
  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
  Let K be the number of rows or columns of a distributed matrix, and
  assume that its process grid has dimension r x c. LOCr( K ) denotes
  the number of elements of K that a process would receive if K were
  distributed over the r processes of its process column. Similarly,
  LOCc( K ) denotes the number of elements of K that a process would
  receive if K were distributed over the c processes of its process
  row. The values of LOCr() and LOCc() may be determined via a call
  to the ScaLAPACK tool function, NUMROC:
          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
  An upper bound for these quantities may be computed by:
          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
  Arguments
  =========
          MP = number of local rows in A and U
          NQ = number of local columns in A and VT
          SIZE = min( M, N )
          SIZEQ = number of local columns in U
          SIZEP = number of local rows in VT
  JOBU    (global input) CHARACTER*1
          Specifies options for computing U:
          = 'V':  the first SIZE columns of U (the left singular
                  vectors) are returned in the array U;
          = 'N':  no columns of U (no left singular vectors) are
                  computed.
  JOBVT   (global input) CHARACTER*1
          Specifies options for computing V**T:
          = 'V':  the first SIZE rows of V**T (the right singular
                  vectors) are returned in the array VT;
          = 'N':  no rows of V**T (no right singular vectors) are
                  computed.
  M       (global input) INTEGER
          The number of rows of the input matrix A.  M >= 0.
  N       (global input) INTEGER
          The number of columns of the input matrix A.  N >= 0.
  A       (local input/workspace) block cyclic DOUBLE PRECISION
          array,
          global dimension (M, N), local dimension (MP, NQ)
          On exit, the contents of A are destroyed.
  IA      (global input) INTEGER
          The row index in the global array A indicating the first
          row of sub( A ).
  JA      (global input) INTEGER
          The column index in the global array A indicating the
          first column of sub( A ).
  DESCA   (global input) INTEGER array of dimension DLEN_
          The array descriptor for the distributed matrix A.
  S       (global output) DOUBLE PRECISION   array, dimension SIZE
          The singular values of A, sorted so that S(i) >= S(i+1).
  U       (local output) DOUBLE PRECISION   array, local dimension
          (MP, SIZEQ), global dimension (M, SIZE)
          if JOBU = 'V', U contains the first min(m,n) columns of U
          if JOBU = 'N', U is not referenced.
  IU      (global input) INTEGER
          The row index in the global array U indicating the first
          row of sub( U ).
  JU      (global input) INTEGER
          The column index in the global array U indicating the
          first column of sub( U ).
  DESCU   (global input) INTEGER array of dimension DLEN_
          The array descriptor for the distributed matrix U.
  VT      (local output) DOUBLE PRECISION   array, local dimension
          (SIZEP, NQ), global dimension (SIZE, N).
          If JOBVT = 'V', VT contains the first SIZE rows of
          V**T. If JOBVT = 'N', VT is not referenced.
  IVT     (global input) INTEGER
          The row index in the global array VT indicating the first
          row of sub( VT ).
  JVT     (global input) INTEGER
          The column index in the global array VT indicating the
          first column of sub( VT ).
  DESCVT   (global input) INTEGER array of dimension DLEN_
          The array descriptor for the distributed matrix VT.
  WORK    (local workspace/output) DOUBLE PRECISION   array, dimension
          (LWORK)
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  LWORK   (local input) INTEGER
          The dimension of the array WORK.
          LWORK >= 1 + 6*SIZEB + MAX(WATOBD, WBDTOSVD),
          where SIZEB = MAX(M,N), and WATOBD and WBDTOSVD refer,
          respectively, to the workspace required to bidiagonalize
          the matrix A and to go from the bidiagonal matrix to the
          singular value decomposition U*S*VT.
          For WATOBD, the following holds:
          WATOBD = MAX(MAX(WPDLANGE,WPDGEBRD),
                       MAX(WPDLARED2D,WP(pre)LARED1D)),
          where WPDLANGE, WPDLARED1D, WPDLARED2D, WPDGEBRD are the
          workspaces required respectively for the subprograms
          PDLANGE, PDLARED1D, PDLARED2D, PDGEBRD. Using the
          standard notation
          MP = NUMROC( M, MB, MYROW, DESCA( CTXT_ ), NPROW),
          NQ = NUMROC( N, NB, MYCOL, DESCA( LLD_ ), NPCOL),
          the workspaces required for the above subprograms are
          WPDLANGE = MP,
          WPDLARED1D = NQ0,
          WPDLARED2D = MP0,
          WPDGEBRD = NB*(MP + NQ + 1) + NQ,
          where NQ0 and MP0 refer, respectively, to the values obtained
          at MYCOL = 0 and MYROW = 0. In general, the upper limit for
          the workspace is given by a workspace required on
          processor (0,0):
          WATOBD <= NB*(MP0 + NQ0 + 1) + NQ0.
          In case of a homogeneous process grid this upper limit can
          be used as an estimate of the minimum workspace for every
          processor.
          For WBDTOSVD, the following holds:
          WBDTOSVD = SIZE*(WANTU*NRU + WANTVT*NCVT) +
                     MAX(WDBDSQR,
                         MAX(WANTU*WPDORMBRQLN, WANTVT*WPDORMBRPRT)),
          where
                          1, if left(right) singular vectors are wanted
          WANTU(WANTVT) =
                          0, otherwise
          and WDBDSQR, WPDORMBRQLN and WPDORMBRPRT refer respectively
          to the workspace required for the subprograms DBDSQR,
          PDORMBR(QLN), and PDORMBR(PRT), where QLN and PRT are the
          values of the arguments VECT, SIDE, and TRANS in the call
          to PDORMBR. NRU is equal to the local number of rows of
          the matrix U when distributed 1-dimensional "column" of
          processes. Analogously, NCVT is equal to the local number
          of columns of the matrix VT when distributed across
          1-dimensional "row" of processes. Calling the LAPACK
          procedure DBDSQR requires
          WDBDSQR = MAX(1, 4*SIZE )
          on every processor. Finally,
          WPDORMBRQLN = MAX( (NB*(NB-1))/2, (SIZEQ+MP)*NB)+NB*NB,
          WPDORMBRPRT = MAX( (MB*(MB-1))/2, (SIZEP+NQ)*MB )+MB*MB,
          If LWORK = -1, then LWORK is global input and a workspace
          query is assumed; the routine only calculates the minimum
          size for the work array. The required workspace is returned
          as the first element of WORK and no error message is issued
          by PXERBLA.
  INFO    (output) INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value

Display dynamic version

FontSize: - +

001        SUBROUTINE PDGESVD(JOBU , JOBVT , M , N , A , IA , JA , DESCA , S , U , IU , JU , DESCU ,
002       +VT , IVT , JVT , DESCVT , WORK , LWORK , INFO)
003  
004  *     -- ScaLAPACK routine(version 1.7) --
005  *     Univ. of Tennessee , Oak Ridge National Laboratory
006  *     and Univ. of California Berkeley.
007  *     Jan 2006
008  
009  *     > 0 : if DBDSQR did not converge
010  *     If INFO = MIN(M , N) + 1 , then PDGESVD has detected
011  *     heterogeneity by finding that eigenvalues were not
012  *     identical across the process grid. In this case , the
013  *     accuracy of the results from PDGESVD cannot be
014  *     guaranteed.
015  
016  *     === ==================================================================
017  
018  *     The results of PDGEBRD , and therefore PDGESVD , may vary slightly
019  *     from run to run with the same input data. If repeatability is an
020  *     issue , call BLACS_SET with the appropriate option after defining
021  *     the process grid.
022  
023  *     Alignment requirements
024  *     === ===================
025  
026  *     The routine PDGESVD inherits the same alignement requirement as
027  *     the routine PDGEBRD , namely :
028  
029  *     The distributed submatrix sub( A ) must verify some alignment proper -
030  *     ties , namely the following expressions should be true :
031  *     ( MB_A.EQ.NB_A .AND. IROFFA.EQ.ICOFFA )
032  *     where NB = MB_A = NB_A ,
033  *     IROFFA = MOD( IA - 1 , NB ) , ICOFFA = MOD( JA - 1 , NB ) ,
034  
035  *     === ==================================================================
036  
037  *     .. Parameters ..
038        INTEGER BLOCK_CYCLIC_2D , DLEN_ , DTYPE_ , CTXT_ , M_ , N_ , MB_ , NB_ , RSRC_ ,
039       +CSRC_ , LLD_ , ITHVAL
040        PARAMETER(BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 , CTXT_ = 2 , M_ = 3 , N_ = 4 ,
041       +MB_ = 5 , NB_ = 6 , RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 , ITHVAL = 10)
042        DOUBLE PRECISION ZERO , ONE
043        PARAMETER(ZERO =(0.0D + 0) , ONE =(1.0D + 0))
044  *     ..
045  *     .. Local Scalars ..
046        CHARACTER UPLO
047        INTEGER CONTEXTC , CONTEXTR , I , INDD , INDD2 , INDE , INDE2 , INDTAUP , INDTAUQ ,
048       +INDU , INDV , INDWORK , IOFFD , IOFFE , ISCALE , J , K , LDU , LDVT , LLWORK ,
049       +LWMIN , MAXIM , MB , MP , MYPCOL , MYPCOLC , MYPCOLR , MYPROW , MYPROWC ,
050       +MYPROWR , NB , NCVT , NPCOL , NPCOLC , NPCOLR , NPROCS , NPROW , NPROWC ,
051       +NPROWR , NQ , NRU , SIZE , SIZEB , SIZEP , SIZEPOS , SIZEQ , WANTU , WANTVT ,
052       +WATOBD , WBDTOSVD , WDBDSQR , WPDGEBRD , WPDLANGE , WPDORMBRPRT ,
053       +WPDORMBRQLN
054        DOUBLE PRECISION ANRM , BIGNUM , EPS , RMAX , RMIN , SAFMIN , SIGMA , SMLNUM
055  *     ..
056  *     .. Local Arrays ..
057        INTEGER DESCTU(DLEN_) , DESCTVT(DLEN_) , IDUM1(3) , IDUM2(3)
058        DOUBLE PRECISION C(1 , 1)
059  *     ..
060  *     .. External Functions ..
061        LOGICAL LSAME
062        INTEGER NUMROC
063        DOUBLE PRECISION PDLAMCH , PDLANGE
064        EXTERNAL LSAME , NUMROC , PDLAMCH , PZLANGE  
065  *     ..
066  *     .. External Subroutines ..
067        EXTERNAL BLACS_GET , BLACS_GRIDEXIT , BLACS_GRIDINFO , BLACS_GRIDINIT ,
068       +CHK1MAT , DBDSQR , DESCINIT , DGAMN2D , DGAMX2D , DSCAL , IGAMX2D ,
069       +IGEBR2D , IGEBS2D , PCHK1MAT , PDGEBRD , PDGEMR2D , PDLARED1D ,
070       +PDLARED2D , PDLASCL , PDLASET , PDORMBR , PXERBLA
071  *     ..
072  *     .. Intrinsic Functions ..
073        INTRINSIC MAX , MIN , SQRT , DBLE
074  *     ..
075  *     .. Executable Statements ..
076  *     This is just to keep ftnchek happy
077        IF(BLOCK_CYCLIC_2D*DTYPE_*LLD_*MB_*M_*NB_*N_.LT.0) RETURN
078  
079        CALL BLACS_GRIDINFO(DESCA(CTXT_) , NPROW , NPCOL , MYPROW , MYPCOL)
080        ISCALE = 0
081        INFO = 0
082  
083        IF(NPROW.EQ. - 1) THEN
083               
084            INFO = - (800 + CTXT_)
085        ELSE
086  
086               
087            SIZE = MIN(M , N)
088            SIZEB = MAX(M , N)
089            NPROCS = NPROW*NPCOL
090            IF(M.GE.N) THEN
090                   
091                IOFFD = JA - 1
092                IOFFE = IA - 1
093                SIZEPOS = 1
094            ELSE
094                   
095                IOFFD = IA - 1
096                IOFFE = JA - 1
097                SIZEPOS = 3
098            END IF
099  
100            IF(LSAME(JOBU , 'V')) THEN
100                   
101                WANTU = 1
102            ELSE
102                   
103                WANTU = 0
104            END IF
105            IF(LSAME(JOBVT , 'V')) THEN
105                   
106                WANTVT = 1
107            ELSE
107                   
108                WANTVT = 0
109            END IF
110  
111            CALL CHK1MAT(M , 3 , N , 4 , IA , JA , DESCA , 8 , INFO)
112            IF(WANTU.EQ.1) THEN
112                   
113                CALL CHK1MAT(M , 3 , SIZE , SIZEPOS , IU , JU , DESCU , 13 , INFO)
114            END IF
115            IF(WANTVT.EQ.1) THEN
115                   
116                CALL CHK1MAT(SIZE , SIZEPOS , N , 4 , IVT , JVT , DESCVT , 17 , INFO)
117            END IF
118            CALL IGAMX2D(DESCA(CTXT_) , 'A' , ' ' , 1 , 1 , INFO , 1 , 1 , 1 ,- 1 ,- 1 , 0)
119  
120            IF(INFO.EQ.0) THEN
121  
122  *             Set up pointers into the WORK array.
123  
123                   
124                INDD = 2
125                INDE = INDD + SIZEB + IOFFD
126                INDD2 = INDE + SIZEB + IOFFE
127                INDE2 = INDD2 + SIZEB + IOFFD
128  
129                INDTAUQ = INDE2 + SIZEB + IOFFE
130                INDTAUP = INDTAUQ + SIZEB + JA - 1
131                INDWORK = INDTAUP + SIZEB + IA - 1
132                LLWORK = LWORK - INDWORK + 1
133  
134  *             Initialize contexts for "column" and "row" process matrices.
135  
136                CALL BLACS_GET(DESCA(CTXT_) , 10 , CONTEXTC)
137                CALL BLACS_GRIDINIT(CONTEXTC , 'R' , NPROCS , 1)
138                CALL BLACS_GRIDINFO(CONTEXTC , NPROWC , NPCOLC , MYPROWC ,
139       +        MYPCOLC)
140                CALL BLACS_GET(DESCA(CTXT_) , 10 , CONTEXTR)
141                CALL BLACS_GRIDINIT(CONTEXTR , 'R' , 1 , NPROCS)
142                CALL BLACS_GRIDINFO(CONTEXTR , NPROWR , NPCOLR , MYPROWR ,
143       +        MYPCOLR)
144  
145  *             Set local dimensions of matrices(this is for MB = NB = 1).
146  
147                NRU = NUMROC(M , 1 , MYPROWC , 0 , NPROCS)
148                NCVT = NUMROC(N , 1 , MYPCOLR , 0 , NPROCS)
149                NB = DESCA(NB_)
150                MB = DESCA(MB_)
151                MP = NUMROC(M , MB , MYPROW , DESCA(RSRC_) , NPROW)
152                NQ = NUMROC(N , NB , MYPCOL , DESCA(CSRC_) , NPCOL)
153                IF(WANTVT.EQ.1) THEN
153                       
154                    SIZEP = NUMROC(SIZE , DESCVT(MB_) , MYPROW , DESCVT(RSRC_) ,
155       +            NPROW)
156                ELSE
156                       
157                    SIZEP = 0
158                END IF
159                IF(WANTU.EQ.1) THEN
159                       
160                    SIZEQ = NUMROC(SIZE , DESCU(NB_) , MYPCOL , DESCU(CSRC_) ,
161       +            NPCOL)
162                ELSE
162                       
163                    SIZEQ = 0
164                END IF
165  
166  *             Transmit MAX(NQ0 , MP0).
167  
168                IF(MYPROW.EQ.0 .AND. MYPCOL.EQ.0) THEN
168                       
169                    MAXIM = MAX(NQ , MP)
170                    CALL IGEBS2D(DESCA(CTXT_) , 'All' , ' ' , 1 , 1 , MAXIM , 1)
171                ELSE
171                       
172                    CALL IGEBR2D(DESCA(CTXT_) , 'All' , ' ' , 1 , 1 , MAXIM , 1 , 0 , 0)
173                END IF
174  
175                WPDLANGE = MP
176                WPDGEBRD = NB* (MP + NQ + 1) + NQ
177                WATOBD = MAX(MAX(WPDLANGE , WPDGEBRD) , MAXIM)
178  
179                WDBDSQR = MAX(1 , 4*SIZE)
180                WPDORMBRQLN = MAX((NB* (NB - 1)) / 2 ,(SIZEQ + MP)*NB) + NB*NB
181                WPDORMBRPRT = MAX((MB* (MB - 1)) / 2 ,(SIZEP + NQ)*MB) + MB*MB
182                WBDTOSVD = SIZE* (WANTU*NRU + WANTVT*NCVT) +
183       +        MAX(WDBDSQR , MAX(WANTU*WPDORMBRQLN ,
184       +        WANTVT*WPDORMBRPRT))
185  
186  *             Finally , calculate required workspace.
187  
188                LWMIN = 1 + 6*SIZEB + MAX(WATOBD , WBDTOSVD)
189                WORK(1) = DBLE(LWMIN)
190  
191                IF(WANTU.NE.1 .AND. .NOT.(LSAME(JOBU , 'N'))) THEN
191                       
192                    INFO = - 1
193                ELSE IF(WANTVT.NE.1 .AND. .NOT.(LSAME(JOBVT , 'N'))) THEN
193                       
194                    INFO = - 2
195                ELSE IF(LWORK.LT.LWMIN .AND. LWORK.NE. - 1) THEN
195                       
196                    INFO = - 19
197                END IF
198  
199            END IF
200  
201            IDUM1(1) = WANTU
202            IDUM1(2) = WANTVT
203            IF(LWORK.EQ. - 1) THEN
203                   
204                IDUM1(3) = - 1
205            ELSE
205                   
206                IDUM1(3) = 1
207            END IF
208            IDUM2(1) = 1
209            IDUM2(2) = 2
210            IDUM2(3) = 19
211            CALL PCHK1MAT(M , 3 , N , 4 , IA , JA , DESCA , 8 , 3 , IDUM1 , IDUM2 , INFO)
212            IF(INFO.EQ.0) THEN
212                   
213                IF(WANTU.EQ.1) THEN
213                       
214                    CALL PCHK1MAT(M , 3 , SIZE , 4 , IU , JU , DESCU , 13 , 0 , IDUM1 , IDUM2 ,
215       +            INFO)
216                END IF
217                IF(WANTVT.EQ.1) THEN
217                       
218                    CALL PCHK1MAT(SIZE , 3 , N , 4 , IVT , JVT , DESCVT , 17 , 0 , IDUM1 ,
219       +            IDUM2 , INFO)
220                END IF
221            END IF
222  
223        END IF
224  
225        IF(INFO.NE.0) THEN
225               
226            CALL PXERBLA(DESCA(CTXT_) , 'PDGESVD' ,- INFO)
227            RETURN
228        ELSE IF(LWORK.EQ. - 1) THEN
228               
229            GO TO 40
230        END IF
231  
232  *     Quick return if possible.
233  
234        IF(M.LE.0 .OR. N.LE.0) GO TO 40
235  
236  *     Get machine constants.
237  
238        SAFMIN = PDLAMCH(DESCA(CTXT_) , 'Safe minimum')
239        EPS = PDLAMCH(DESCA(CTXT_) , 'Precision')
240        SMLNUM = SAFMIN / EPS
241        BIGNUM = ONE / SMLNUM
242        RMIN = SQRT(SMLNUM)
243        RMAX = MIN(SQRT(BIGNUM) , ONE / SQRT(SQRT(SAFMIN)))
244  
245  *     Scale matrix to allowable range , if necessary.
246  
247        ANRM = PDLANGE('1' , M , N , A , IA , JA , DESCA , WORK(INDWORK))
248        IF(ANRM.GT.ZERO .AND. ANRM.LT.RMIN) THEN
248               
249            ISCALE = 1
250            SIGMA = RMIN / ANRM
251        ELSE IF(ANRM.GT.RMAX) THEN
251               
252            ISCALE = 1
253            SIGMA = RMAX / ANRM
254        END IF
255  
256        IF(ISCALE.EQ.1) THEN
256               
257            CALL PDLASCL ('G' , ONE , SIGMA , M , N , A , IA , JA , DESCA , INFO)
258        END IF
259  
260        CALL PDGEBRD (M , N , A , IA , JA , DESCA , WORK(INDD) , WORK(INDE) ,
261       +WORK(INDTAUQ) , WORK(INDTAUP) , WORK(INDWORK) , LLWORK ,
262       +INFO)
263  
264  *     Copy D and E to all processes.
265  *     Array D is in local array of dimension :
266  *     LOCc(JA + MIN(M , N) - 1) if M >= N ; LOCr(IA + MIN(M , N) - 1) otherwise.
267  *     Array E is in local array of dimension
268  *     LOCr(IA + MIN(M , N) - 1) if M >= N ; LOCc(JA + MIN(M , N) - 2) otherwise.
269  
270        IF(M.GE.N) THEN
271  *         Distribute D
271               
272            CALL PDLARED1D (N + IOFFD , IA , JA , DESCA , WORK(INDD) , WORK(INDD2) ,
273       +    WORK(INDWORK) , LLWORK)
274  *         Distribute E
275            CALL PDLARED2D (M + IOFFE , IA , JA , DESCA , WORK(INDE) , WORK(INDE2) ,
276       +    WORK(INDWORK) , LLWORK)
277        ELSE
278  *         Distribute D
278               
279            CALL PDLARED2D (M + IOFFD , IA , JA , DESCA , WORK(INDD) , WORK(INDD2) ,
280       +    WORK(INDWORK) , LLWORK)
281  *         Distribute E
282            CALL PDLARED1D (N + IOFFE , IA , JA , DESCA , WORK(INDE) , WORK(INDE2) ,
283       +    WORK(INDWORK) , LLWORK)
284        END IF
285  
286  *     Prepare for calling PDBDSQR.
287  
288        IF(M.GE.N) THEN
288               
289            UPLO = 'U'
290        ELSE
290               
291            UPLO = 'L'
292        END IF
293  
294        INDU = INDWORK
295        INDV = INDU + SIZE*NRU*WANTU
296        INDWORK = INDV + SIZE*NCVT*WANTVT
297  
298        LDU = MAX(1 , NRU)
299        LDVT = MAX(1 , SIZE)
300  
301        CALL DESCINIT(DESCTU , M , SIZE , 1 , 1 , 0 , 0 , CONTEXTC , LDU , INFO)
302        CALL DESCINIT(DESCTVT , SIZE , N , 1 , 1 , 0 , 0 , CONTEXTR , LDVT , INFO)
303  
304        IF(WANTU.EQ.1) THEN
304               
305            CALL PDLASET ('Full' , M , SIZE , ZERO , ONE , WORK(INDU) , 1 , 1 , DESCTU)
306        ELSE
306               
307            NRU = 0
308        END IF
309  
310        IF(WANTVT.EQ.1) THEN
310               
311            CALL PDLASET ('Full' , SIZE , N , ZERO , ONE , WORK(INDV) , 1 , 1 , DESCTVT)
312        ELSE
312               
313            NCVT = 0
314        END IF
315  
316        CALL DBDSQR(UPLO , SIZE , NCVT , NRU , 0 , WORK(INDD2 + IOFFD) ,
317       +WORK(INDE2 + IOFFE) , WORK(INDV) , SIZE , WORK(INDU) , LDU , C , 1 ,
318       +WORK(INDWORK) , INFO)
319  
320  *     Redistribute elements of U and VT in the block - cyclic fashion.
321  
322        IF(WANTU.EQ.1) CALL PDGEMR2D(M , SIZE , WORK(INDU) , 1 , 1 , DESCTU , U , IU ,
323       +JU , DESCU , DESCU(CTXT_))
324  
325        IF(WANTVT.EQ.1) CALL PDGEMR2D(SIZE , N , WORK(INDV) , 1 , 1 , DESCTVT , VT ,
326       +IVT , JVT , DESCVT , DESCVT(CTXT_))
327  
328  *     Set to ZERO "non-square" elements of the larger matrices U , VT.
329  
330        IF(M.GT.N .AND. WANTU.EQ.1) THEN
330               
331            CALL PDLASET ('Full' , M - SIZE , SIZE , ZERO , ZERO , U , IA + SIZE , JU , DESCU)
332        ELSE IF(N.GT.M .AND. WANTVT.EQ.1) THEN
332               
333            CALL PDLASET ('Full' , SIZE , N - SIZE , ZERO , ZERO , VT , IVT , JVT + SIZE ,
334       +    DESCVT)
335        END IF
336  
337  *     Multiply Householder rotations from bidiagonalized matrix.
338  
339        IF(WANTU.EQ.1) CALL PDORMBR ('Q' , 'L' , 'N' , M , SIZE , N , A , IA , JA , DESCA ,
340       +WORK(INDTAUQ) , U , IU , JU , DESCU ,
341       +WORK(INDWORK) , LLWORK , INFO)
342  
343        IF(WANTVT.EQ.1) CALL PDORMBR ('P' , 'R' , 'T' , SIZE , N , M , A , IA , JA , DESCA ,
344       +WORK(INDTAUP) , VT , IVT , JVT , DESCVT ,
345       +WORK(INDWORK) , LLWORK , INFO)
346  
347  *     Copy singular values into output array S.
348  
349        DO 10 I = 1 , SIZE
349               
350            S(I) = WORK(INDD2 + IOFFD + I - 1)
351     10 CONTINUE
352  
353  *     If matrix was scaled , then rescale singular values appropriately.
354  
355        IF(ISCALE.EQ.1) THEN
355               
356            CALL DSCAL(SIZE , ONE / SIGMA , S , 1)
357        END IF
358  
359  *     Compare every ith eigenvalue , or all if there are only a few ,
360  *     across the process grid to check for heterogeneity.
361  
362        IF(SIZE.LE.ITHVAL) THEN
362               
363            J = SIZE
364            K = 1
365        ELSE
365               
366            J = SIZE / ITHVAL
367            K = ITHVAL
368        END IF
369  
370        DO 20 I = 1 , J
370               
371            WORK(I + INDE) = S((I - 1)*K + 1)
372            WORK(I + INDD2) = S((I - 1)*K + 1)
373     20 CONTINUE
374  
375        CALL DGAMN2D(DESCA(CTXT_) , 'a' , ' ' , J , 1 , WORK(1 + INDE) , J , 1 , 1 ,- 1 ,- 1 , 0)
376        CALL DGAMX2D(DESCA(CTXT_) , 'a' , ' ' , J , 1 , WORK(1 + INDD2) , J , 1 , 1 ,- 1 ,- 1 , 0)
377  
378        DO 30 I = 1 , J
378               
379            IF((WORK(I + INDE) - WORK(I + INDD2)).NE.ZERO) THEN
379                   
380                INFO = SIZE + 1
381            END IF
382     30 CONTINUE
383  
384     40 CONTINUE
385  
386        CALL BLACS_GRIDEXIT(CONTEXTC)
387        CALL BLACS_GRIDEXIT(CONTEXTR)
388  
389  *     End of PDGESVD
390  
391        RETURN
392        END6572

Variables in Routine PDGESVD()

Summary Report
Data Type	Quantity	Size(byte)
CHARACTER	1	1
DOUBLE PRECISION	13	52
INTEGER	73	312
LOGICAL	1	1
REAL	2	8
TOTAL	90	374

List of Variables

CHARACTER

UPLO

DOUBLE PRECISION

ANRM BIGNUM C(1,1) EPS ONE
PDLAMCH PDLANGE RMAX RMIN SAFMIN
SIGMA SMLNUM ZERO

INTEGER

BLOCK_CYCLIC_2D CONTEXTC CONTEXTR CSRC_ CTXT_
DESCTU(DLEN_) DESCTVT(DLEN_) DLEN_ DTYPE_ I
IDUM1(3) IDUM2(3) INDD INDD2 INDE
INDE2 INDTAUP INDTAUQ INDU INDV
INDWORK INFO IOFFD IOFFE ISCALE
ITHVAL J K LDU LDVT
LLD_ LLWORK LWMIN M_ MAXIM
MB MB_ MP MYPCOL MYPCOLC
MYPCOLR MYPROW MYPROWC MYPROWR N_
NB NB_ NCVT NPCOL NPCOLC
NPCOLR NPROCS NPROW NPROWC NPROWR
NQ NRU NUMROC RSRC_ SIZE
SIZEB SIZEP SIZEPOS SIZEQ WANTU
WANTVT WATOBD WBDTOSVD WDBDSQR WPDGEBRD
WPDLANGE WPDORMBRPRT WPDORMBRQLN

LOGICAL

LSAME

REAL

S WORK

Variables Dependence Graph
Put the mouse over a right hand side variable to display the corresponding line of the dependence
-
- -
ANRM <--- INDWORK, PDLANGE, WORK

BIGNUM <--- ONEBIGNUM = ONE/SMLNUM, SMLNUMBIGNUM = ONE/SMLNUM

EPS <--- PDLAMCHEPS = PDLAMCH(DESCA(CTXT_),'Precision'), CTXT_EPS = PDLAMCH(DESCA(CTXT_),'Precision')

I <--- JDO 20 I = 1,J{2DO 30 I = 1,J}, SIZEDO 10 I = 1,SIZE

IDUM1 <--- WANTUIDUM1(1) = WANTU, WANTVTIDUM1(2) = WANTVT

INDD2 <--- INDEINDD2 = INDE + SIZEB + IOFFE, IOFFEINDD2 = INDE + SIZEB + IOFFE, SIZEBINDD2 = INDE + SIZEB + IOFFE

INDE <--- INDDINDE = INDD + SIZEB + IOFFD, IOFFDINDE = INDD + SIZEB + IOFFD, SIZEBINDE = INDD + SIZEB + IOFFD

INDE2 <--- INDD2INDE2 = INDD2 + SIZEB + IOFFD, IOFFDINDE2 = INDD2 + SIZEB + IOFFD, SIZEBINDE2 = INDD2 + SIZEB + IOFFD

INDTAUP <--- INDTAUQINDTAUP = INDTAUQ + SIZEB + JA - 1, SIZEBINDTAUP = INDTAUQ + SIZEB + JA - 1

INDTAUQ <--- INDE2INDTAUQ = INDE2 + SIZEB + IOFFE, IOFFEINDTAUQ = INDE2 + SIZEB + IOFFE, SIZEBINDTAUQ = INDE2 + SIZEB + IOFFE

INDU <--- INDWORKINDU = INDWORK

INDV <--- INDUINDV = INDU + SIZE*NRU*WANTU, NRUINDV = INDU + SIZE*NRU*WANTU, SIZEINDV = INDU + SIZE*NRU*WANTU, WANTUINDV = INDU + SIZE*NRU*WANTU

INDWORK <--- INDTAUPINDWORK = INDTAUP + SIZEB + IA - 1, INDVINDWORK = INDV + SIZE*NCVT*WANTVT, NCVTINDWORK = INDV + SIZE*NCVT*WANTVT, SIZEINDWORK = INDV + SIZE*NCVT*WANTVT, SIZEBINDWORK = INDTAUP + SIZEB + IA - 1, WANTVTINDWORK = INDV + SIZE*NCVT*WANTVT

INFO <--- SIZEINFO = SIZE + 1, CTXT_INFO = - (800+CTXT_)

J <--- ITHVALJ = SIZE/ITHVAL, SIZEJ = SIZE{2J = SIZE/ITHVAL}

K <--- ITHVALK = ITHVAL

LDU <--- NRULDU = MAX(1,NRU)

LDVT <--- SIZELDVT = MAX(1,SIZE)

LLWORK <--- INDWORKLLWORK = LWORK - INDWORK + 1

LWMIN <--- SIZEB, WATOBD, WBDTOSVD

MAXIM <--- MPMAXIM = MAX(NQ,MP), NQMAXIM = MAX(NQ,MP)

MB <--- MB_MB = DESCA(MB_)

MP <--- MB, MYPROW, NPROW, NUMROC, RSRC_

NB <--- NB_NB = DESCA(NB_)

NCVT <--- MYPCOLRNCVT = NUMROC(N,1,MYPCOLR,0,NPROCS), NPROCSNCVT = NUMROC(N,1,MYPCOLR,0,NPROCS), NUMROCNCVT = NUMROC(N,1,MYPCOLR,0,NPROCS)

NPROCS <--- NPCOLNPROCS = NPROW*NPCOL, NPROWNPROCS = NPROW*NPCOL

NQ <--- MYPCOL, NB, NPCOL, NUMROC, CSRC_

NRU <--- MYPROWCNRU = NUMROC(M,1,MYPROWC,0,NPROCS), NPROCSNRU = NUMROC(M,1,MYPROWC,0,NPROCS), NUMROCNRU = NUMROC(M,1,MYPROWC,0,NPROCS)

RMAX <--- BIGNUM, ONE, SAFMIN

RMIN <--- SMLNUMRMIN = SQRT(SMLNUM)

S <--- IS(I) = WORK(INDD2+IOFFD+I-1), INDD2S(I) = WORK(INDD2+IOFFD+I-1), IOFFDS(I) = WORK(INDD2+IOFFD+I-1), WORKS(I) = WORK(INDD2+IOFFD+I-1)

SAFMIN <--- PDLAMCH, CTXT_

SIGMA <--- ANRMSIGMA = RMIN/ANRM{2SIGMA = RMAX/ANRM}, RMAXSIGMA = RMAX/ANRM, RMINSIGMA = RMIN/ANRM

SIZEP <--- MB_, MYPROW, NPROW, NUMROC, RSRC_, SIZE

SIZEQ <--- MYPCOL, NB_, NPCOL, NUMROC, CSRC_, SIZE

SMLNUM <--- EPSSMLNUM = SAFMIN/EPS, SAFMINSMLNUM = SAFMIN/EPS

WATOBD <--- MAXIM, WPDGEBRD, WPDLANGE

WBDTOSVD <--- NCVT, NRU, SIZE, WANTU, WANTVT, WDBDSQR, WPDORMBRPRT, WPDORMBRQLN

WDBDSQR <--- SIZEWDBDSQR = MAX(1,4*SIZE)

WORK <--- IWORK(I+INDE) = S((I-1)*K+1){2WORK(I+INDD2) = S((I-1)*K+1)}, KWORK(I+INDE) = S((I-1)*K+1){2WORK(I+INDD2) = S((I-1)*K+1)}, LWMINWORK(1) = DBLE(LWMIN), SWORK(I+INDE) = S((I-1)*K+1){2WORK(I+INDD2) = S((I-1)*K+1)}

WPDGEBRD <--- MPWPDGEBRD = NB* (MP+NQ+1) + NQ, NBWPDGEBRD = NB* (MP+NQ+1) + NQ, NQWPDGEBRD = NB* (MP+NQ+1) + NQ

WPDLANGE <--- MPWPDLANGE = MP

WPDORMBRPRT <--- MB, NQ, SIZEP

WPDORMBRQLN <--- MP, NB, SIZEQ