Routine: PCHEEVD()  File: SRC\pcheevd.f

 
 
# lines: 441
  # code: 441
  # comment: 0
  # blank:0
# Variables:93
# Callers:0
# Callings:9
# Words:257
# Keywords:158
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PCHEEVD computes all the eigenvalues and eigenvectors of a Hermitian
  matrix A by using a divide and conquer algorithm.
  Arguments
  =========
     NP = the number of rows local to a given process.
     NQ = the number of columns local to a given process.
  JOBZ    (input) CHARACTER*1
          = 'N':  Compute eigenvalues only;    (NOT IMPLEMENTED YET)
          = 'V':  Compute eigenvalues and eigenvectors.
  UPLO    (global input) CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored:
          = 'U':  Upper triangular
          = 'L':  Lower triangular
  N       (global input) INTEGER
          The number of rows and columns of the matrix A.  N >= 0.
  A       (local input/workspace) block cyclic COMPLEX array,
          global dimension (N, N), local dimension ( LLD_A,
          LOCc(JA+N-1) )
          On entry, the symmetric matrix A.  If UPLO = 'U', only the
          upper triangular part of A is used to define the elements of
          the symmetric matrix.  If UPLO = 'L', only the lower
          triangular part of A is used to define the elements of the
          symmetric matrix.
          On exit, the lower triangle (if UPLO='L') or the upper
          triangle (if UPLO='U') of A, including the diagonal, is
          destroyed.
  IA      (global input) INTEGER
          A's global row index, which points to the beginning of the
          submatrix which is to be operated on.
  JA      (global input) INTEGER
          A's global column index, which points to the beginning of
          the submatrix which is to be operated on.
  DESCA   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix A.
          If DESCA( CTXT_ ) is incorrect, PCHEEVD cannot guarantee
          correct error reporting.
  W       (global output) REAL array, dimension (N)
          If INFO=0, the eigenvalues in ascending order.
  Z       (local output) COMPLEX array,
          global dimension (N, N),
          local dimension ( LLD_Z, LOCc(JZ+N-1) )
          Z contains the orthonormal eigenvectors of the matrix A.
  IZ      (global input) INTEGER
          Z's global row index, which points to the beginning of the
          submatrix which is to be operated on.
  JZ      (global input) INTEGER
          Z's global column index, which points to the beginning of
          the submatrix which is to be operated on.
  DESCZ   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix Z.
          DESCZ( CTXT_ ) must equal DESCA( CTXT_ )
  WORK    (local workspace/output) COMPLEX array,
          dimension (LWORK)
          On output, WORK(1) returns the workspace needed for the
          computation.
  LWORK   (local input) INTEGER
          If eigenvectors are requested:
            LWORK = N + ( NP0 + MQ0 + NB ) * NB,
          with  NP0 = NUMROC( MAX( N, NB, 2 ), NB, 0, 0, NPROW )
                MQ0 = NUMROC( MAX( N, NB, 2 ), NB, 0, 0, NPCOL )
          If LWORK = -1, then LWORK is global input and a workspace
          query is assumed; the routine calculates the size for all
          work arrays. Each of these values is returned in the first
          entry of the corresponding work array, and no error message
          is issued by PXERBLA.
  RWORK   (local workspace/output) REAL array,
          dimension (LRWORK)
          On output RWORK(1) returns the real workspace needed to
          guarantee completion.  If the input parameters are incorrect,
          RWORK(1) may also be incorrect.
  LRWORK  (local input) INTEGER
          Size of RWORK array.
          LRWORK >= 1 + 9*N + 3*NP*NQ,
          NP = NUMROC( N, NB, MYROW, IAROW, NPROW )
          NQ = NUMROC( N, NB, MYCOL, IACOL, NPCOL )
  IWORK   (local workspace/output) INTEGER array, dimension (LIWORK)
          On output IWORK(1) returns the integer workspace needed.
  LIWORK  (input) INTEGER
          The dimension of the array IWORK.
          LIWORK = 7*N + 8*NPCOL + 2
  INFO    (global output) INTEGER
          = 0:  successful exit
          < 0:  If the i-th argument is an array and the j-entry had
                an illegal value, then INFO = -(i*100+j), if the i-th
                argument is a scalar and had an illegal value, then
                INFO = -i.
          > 0:  If INFO = 1 through N, the i(th) eigenvalue did not
                converge in PSLAED3.
  Alignment requirements
  ======================
  The distributed submatrices sub( A ), sub( Z ) must verify
  some alignment properties, namely the following expression
  should be true:
  ( MB_A.EQ.NB_A.EQ.MB_Z.EQ.NB_Z .AND. IROFFA.EQ.ICOFFA .AND.
    IROFFA.EQ.0 .AND.IROFFA.EQ.IROFFZ. AND. IAROW.EQ.IZROW)
    with IROFFA = MOD( IA-1, MB_A )
     and ICOFFA = MOD( JA-1, NB_A ).
  Further Details
  ======= =======
  Contributed by Francoise Tisseur, University of Manchester.
  Reference:  F. Tisseur and J. Dongarra, "A Parallel Divide and
              Conquer Algorithm for the Symmetric Eigenvalue Problem
              on Distributed Memory Architectures",
              SIAM J. Sci. Comput., 6:20 (1999), pp. 2223--2236.
              (see also LAPACK Working Note 132)
                http://www.netlib.org/lapack/lawns/lawn132.ps
  =====================================================================
     .. Parameters ..

 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

 
001        SUBROUTINE PCHEEVD( JOBZ , UPLO , N , A , IA , JA , DESCA , W , Z , IZ , JZ ,
002       $DESCZ , WORK , LWORK , RWORK , LRWORK , IWORK ,
003       $LIWORK , INFO )
004  
005  *     -- ScaLAPACK routine(version 1.7) --
006  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
007  *     and University of California , Berkeley.
008  *     March 25 , 2002
009  
010  *     .. Scalar Arguments ..
011        CHARACTER JOBZ , UPLO
012        INTEGER IA , INFO , IZ , JA , JZ , LIWORK , LRWORK , LWORK , N
013        INTEGER BLOCK_CYCLIC_2D , DLEN_ , DTYPE_ , CTXT_ , M_ , N_ ,
014       $MB_ , NB_ , RSRC_ , CSRC_ , LLD_
015        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
016       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
017       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
018        REAL ZERO , ONE
019        PARAMETER( ZERO = 0.0E + 0 , ONE = 1.0E + 0 )
020        COMPLEX CZERO , CONE
021        PARAMETER( CZERO =( 0.0E + 0 , 0.0E + 0 ) ,
022       $CONE =( 1.0E + 0 , 0.0E + 0 ) )
023  *     ..
024  *     .. Local Scalars ..
025        LOGICAL LOWER , LQUERY
026        INTEGER CSRC_A , I , IACOL , IAROW , ICOFFA , IINFO , IIZ ,
027       $INDD , INDE , INDE2 , INDRWORK , INDTAU , INDWORK ,
028       $INDZ , IPR , IPZ , IROFFA , IROFFZ , ISCALE , IZCOL ,
029       $IZROW , J , JJZ , LDR , LDZ , LIWMIN , LLRWORK ,
030       $LLWORK , LRWMIN , LWMIN , MB_A , MYCOL , MYROW , NB ,
031       $NB_A , NN , NP0 , NPCOL , NPROW , NQ , NQ0 , OFFSET ,
032       $RSRC_A
033        REAL ANRM , BIGNUM , EPS , RMAX , RMIN , SAFMIN , SIGMA ,
034       $SMLNUM
035  *     ..
036  *     .. Local Arrays ..
037        INTEGER DESCRZ( 9 ) , IDUM1( 2 ) , IDUM2( 2 )
038  *     ..
039  *     .. External Functions ..
040        LOGICAL LSAME
041        INTEGER INDXG2L , INDXG2P , NUMROC
042        REAL PCLANHE , PSLAMCH
043        EXTERNAL LSAME , INDXG2L , INDXG2P , NUMROC , PCLANHE ,
044       $PSLAMCH
045  *     ..
046  *     .. External Subroutines ..
047        EXTERNAL BLACS_GRIDINFO , CHK1MAT , DESCINIT , INFOG2L ,
048       $PCELGET , PCHETRD , PCHK2MAT , PCLASCL , PCLASET ,
049       $PCUNMTR , PSLARED1D , PSLASET , PSSTEDC , PXERBLA ,
050       $SSCAL
051  *     ..
052  *     .. Intrinsic Functions ..
053        INTRINSIC CMPLX , ICHAR , MAX , MIN , MOD , REAL , SQRT
054  *     ..
055  *     .. Executable Statements ..
056  *     This is just to keep ftnchek and toolpack / 1 happy
057        IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
058       $    RSRC_.LT.0 )RETURN
059  
060            INFO = 0
061  
062  *         Quick return
063  
064            IF( N.EQ.0 )
065       $        RETURN
066  
067  *             Test the input arguments.
068  
069                CALL BLACS_GRIDINFO( DESCA( CTXT_ ) , NPROW , NPCOL , MYROW , MYCOL )
070  
071                IF( NPROW.EQ. - 1 ) THEN
072                    INFO = - ( 700 + CTXT_ )
073                ELSE
074                    CALL CHK1MAT( N , 3 , N , 3 , IA , JA , DESCA , 7 , INFO )
075                    CALL CHK1MAT( N , 3 , N , 3 , IZ , JZ , DESCZ , 12 , INFO )
076                    IF( INFO.EQ.0 ) THEN
077                        LOWER = LSAME( UPLO , 'L' )
078                        NB_A = DESCA( NB_ )
079                        MB_A = DESCA( MB_ )
080                        NB = NB_A
081                        RSRC_A = DESCA( RSRC_ )
082                        CSRC_A = DESCA( CSRC_ )
083                        IROFFA = MOD( IA - 1 , MB_A )
084                        ICOFFA = MOD( JA - 1 , NB_A )
085                        IAROW = INDXG2P( IA , NB_A , MYROW , RSRC_A , NPROW )
086                        IACOL = INDXG2P( JA , MB_A , MYCOL , CSRC_A , NPCOL )
087                        NP0 = NUMROC( N , NB , MYROW , IAROW , NPROW )
088                        NQ0 = NUMROC( N , NB , MYCOL , IACOL , NPCOL )
089                        IROFFZ = MOD( IZ - 1 , MB_A )
090                        CALL INFOG2L( IZ , JZ , DESCZ , NPROW , NPCOL , MYROW , MYCOL ,
091       $                IIZ , JJZ , IZROW , IZCOL )
092                        LQUERY =( LWORK.EQ. - 1 .OR. LIWORK.EQ. - 1 .OR. LRWORK.EQ. - 1 )
093  
094  *                     Compute the total amount of space needed
095  
096                        NN = MAX( N , NB , 2 )
097                        NQ = NUMROC( NN , NB , 0 , 0 , NPCOL )
098                        LWMIN = N + ( NP0 + NQ + NB )*NB
099                        LRWMIN = 1 + 9*N + 3*NP0*NQ0
100                        LIWMIN = 7*N + 8*NPCOL + 2
101                        WORK( 1 ) = CMPLX( LWMIN )
102                        RWORK( 1 ) = REAL( LRWMIN )
103                        IWORK( 1 ) = LIWMIN
104                        IF( .NOT.LSAME( JOBZ , 'V' ) ) THEN
105                            INFO = - 1
106                        ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO , 'U' ) ) ) THEN
107                            INFO = - 2
108                        ELSE IF( LWORK.LT.LWMIN .AND. LWORK.NE. - 1 ) THEN
109                            INFO = - 14
110                        ELSE IF( LRWORK.LT.LRWMIN .AND. LRWORK.NE. - 1 ) THEN
111                            INFO = - 16
112                        ELSE IF( IROFFA.NE.0 ) THEN
113                            INFO = - 4
114                        ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
115                            INFO = - ( 700 + NB_ )
116                        ELSE IF( IROFFA.NE.IROFFZ ) THEN
117                            INFO = - 10
118                        ELSE IF( IAROW.NE.IZROW ) THEN
119                            INFO = - 10
120                        ELSE IF( DESCA( M_ ).NE.DESCZ( M_ ) ) THEN
121                            INFO = - ( 1200 + M_ )
122                        ELSE IF( DESCA( N_ ).NE.DESCZ( N_ ) ) THEN
123                            INFO = - ( 1200 + N_ )
124                        ELSE IF( DESCA( MB_ ).NE.DESCZ( MB_ ) ) THEN
125                            INFO = - ( 1200 + MB_ )
126                        ELSE IF( DESCA( NB_ ).NE.DESCZ( NB_ ) ) THEN
127                            INFO = - ( 1200 + NB_ )
128                        ELSE IF( DESCA( RSRC_ ).NE.DESCZ( RSRC_ ) ) THEN
129                            INFO = - ( 1200 + RSRC_ )
130                        ELSE IF( DESCA( CTXT_ ).NE.DESCZ( CTXT_ ) ) THEN
131                            INFO = - ( 1200 + CTXT_ )
132                        END IF
133                    END IF
134                    IF( LOWER ) THEN
135                        IDUM1( 1 ) = ICHAR( 'L' )
136                    ELSE
137                        IDUM1( 1 ) = ICHAR( 'U' )
138                    END IF
139                    IDUM2( 1 ) = 2
140                    IF( LWORK.EQ. - 1 ) THEN
141                        IDUM1( 2 ) = - 1
142                    ELSE
143                        IDUM1( 2 ) = 1
144                    END IF
145                    IDUM2( 2 ) = 14
146                    CALL PCHK2MAT( N , 3 , N , 3 , IA , JA , DESCA , 7 , N , 3 , N , 3 , IZ ,
147       $            JZ , DESCZ , 12 , 2 , IDUM1 , IDUM2 , INFO )
148                END IF
149  
150                IF( INFO.NE.0 ) THEN
151                    CALL PXERBLA( DESCA( CTXT_ ) , 'PCHEEVD' , - INFO )
152                    RETURN
153                ELSE IF( LQUERY ) THEN
154                    RETURN
155                END IF
156  
157  *             Get machine constants.
158  
159                SAFMIN = PSLAMCH( DESCA( CTXT_ ) , 'Safe minimum' )
160                EPS = PSLAMCH( DESCA( CTXT_ ) , 'Precision' )
161                SMLNUM = SAFMIN / EPS
162                BIGNUM = ONE / SMLNUM
163                RMIN = SQRT( SMLNUM )
164                RMAX = MIN( SQRT( BIGNUM ) , ONE / SQRT( SQRT( SAFMIN ) ) )
165  
166  *             Set up pointers into the WORK array
167  
168                INDTAU = 1
169                INDWORK = INDTAU + N
170                LLWORK = LWORK - INDWORK + 1
171  
172  *             Set up pointers into the RWORK array
173  
174                INDE = 1
175                INDD = INDE + N
176                INDE2 = INDD + N
177                INDRWORK = INDE2 + N
178                LLRWORK = LRWORK - INDRWORK + 1
179  
180  *             Scale matrix to allowable range , if necessary.
181  
182                ISCALE = 0
183  
184                ANRM = PCLANHE( 'M' , UPLO , N , A , IA , JA , DESCA ,
185       $        RWORK( INDRWORK ) )
186  
187                IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
188                    ISCALE = 1
189                    SIGMA = RMIN / ANRM
190                ELSE IF( ANRM.GT.RMAX ) THEN
191                    ISCALE = 1
192                    SIGMA = RMAX / ANRM
193                END IF
194  
195                IF( ISCALE.EQ.1 ) THEN
196                    CALL PCLASCL ( UPLO , ONE , SIGMA , N , N , A , IA , JA , DESCA , IINFO )
197                END IF
198  
199  *             Reduce Hermitian matrix to tridiagonal form.
200  
201                CALL PCHETRD ( UPLO , N , A , IA , JA , DESCA , RWORK( INDD ) ,
202       $        RWORK( INDE2 ) , WORK( INDTAU ) , WORK( INDWORK ) ,
203       $        LLWORK , IINFO )
204  
205  *             Copy the values of D , E to all processes
206  
207  *             Here PxLARED1D is used to redistribute the tridiagonal matrix.
208  *             PxLARED1D , however , doesn't yet work with arbritary matrix
209  *             distributions so we have PxELGET as a backup.
210  
211                OFFSET = 0
212                IF( IA.EQ.1 .AND. JA.EQ.1 .AND. RSRC_A.EQ.0 .AND. CSRC_A.EQ.0 )
213       $            THEN
214                    CALL PSLARED1D ( N , IA , JA , DESCA , RWORK( INDD ) , W ,
215       $            RWORK( INDRWORK ) , LLRWORK )
216  
217                    CALL PSLARED1D ( N , IA , JA , DESCA , RWORK( INDE2 ) ,
218       $            RWORK( INDE ) , RWORK( INDRWORK ) , LLRWORK )
219                    IF( .NOT.LOWER )
220       $                OFFSET = 1
221                    ELSE
222                        DO 10 I = 1 , N
223                            CALL PCELGET( 'A' , ' ' , WORK( INDWORK ) , A , I + IA - 1 , I + JA - 1 ,
224       $                    DESCA )
225                            W( I ) = REAL( WORK( INDWORK ) )
226     10                 CONTINUE
227                        IF( LSAME( UPLO , 'U' ) ) THEN
228                            DO 20 I = 1 , N - 1
229                                CALL PCELGET( 'A' , ' ' , WORK( INDWORK ) , A , I + IA - 1 , I + JA ,
230       $                        DESCA )
231                                RWORK( INDE + I - 1 ) = REAL( WORK( INDWORK ) )
232     20                     CONTINUE
233                        ELSE
234                            DO 30 I = 1 , N - 1
235                                CALL PCELGET( 'A' , ' ' , WORK( INDWORK ) , A , I + IA , I + JA - 1 ,
236       $                        DESCA )
237                                RWORK( INDE + I - 1 ) = REAL( WORK( INDWORK ) )
238     30                     CONTINUE
239                        END IF
240                    END IF
241  
242  *                 Call PSSTEDC to compute eigenvalues and eigenvectors.
243  
244                    INDZ = INDE + N
245                    INDRWORK = INDZ + NP0*NQ0
246                    LLRWORK = LRWORK - INDRWORK + 1
247                    LDR = MAX( 1 , NP0 )
248                    CALL DESCINIT( DESCRZ , DESCZ( M_ ) , DESCZ( N_ ) , DESCZ( MB_ ) ,
249       $            DESCZ( NB_ ) , DESCZ( RSRC_ ) , DESCZ( CSRC_ ) ,
250       $            DESCZ( CTXT_ ) , LDR , INFO )
251                    CALL PCLASET ( 'Full' , N , N , CZERO , CONE , Z , IZ , JZ , DESCZ )
252                    CALL PSLASET ( 'Full' , N , N , ZERO , ONE , RWORK( INDZ ) , 1 , 1 ,
253       $            DESCRZ )
254                    CALL PSSTEDC ( 'I' , N , W , RWORK( INDE + OFFSET ) , RWORK( INDZ ) , IZ ,
255       $            JZ , DESCRZ , RWORK( INDRWORK ) , LLRWORK , IWORK ,
256       $            LIWORK , IINFO )
257  
258                    LDZ = DESCZ( LLD_ )
259                    LDR = DESCRZ( LLD_ )
260                    IIZ = INDXG2L( IZ , NB , MYROW , MYROW , NPROW )
261                    JJZ = INDXG2L( JZ , NB , MYCOL , MYCOL , NPCOL )
262                    IPZ = IIZ + ( JJZ - 1 )*LDZ
263                    IPR = INDZ - 1 + IIZ + ( JJZ - 1 )*LDR
264                    DO 50 J = 0 , NQ0 - 1
265                        DO 40 I = 0 , NP0 - 1
266                            Z( IPZ + I + J*LDZ ) = RWORK( IPR + I + J*LDR )
267     40                 CONTINUE
268     50             CONTINUE
269  
270  *                 Z = Q * Z
271  
272                    CALL PCUNMTR ( 'L' , UPLO , 'N' , N , N , A , IA , JA , DESCA ,
273       $            WORK( INDTAU ) , Z , IZ , JZ , DESCZ , WORK( INDWORK ) ,
274       $            LLWORK , IINFO )
275  
276  *                 If matrix was scaled , then rescale eigenvalues appropriately.
277  
278                    IF( ISCALE.EQ.1 ) THEN
279                        CALL SSCAL( N , ONE / SIGMA , W , 1 )
280                    END IF
281  
282                    WORK( 1 ) = CMPLX( LWMIN )
283                    RWORK( 1 ) = REAL( LRWMIN )
284                    IWORK( 1 ) = LIWMIN
285  
286                    RETURN
287  
288  *                 End of PCHEEVD
289  
290                END