Routine PDPOTF2() Documented By Universal Report

Routine: PDPOTF2()

# lines:	351
# code:	351
# comment:	0
# blank:	0
# Variables:	41
# Callers:	1
# Callings:	0
# Words:	179
# Keywords:	116

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PDPOTF2 computes the Cholesky factorization of a real symmetric
  positive definite distributed matrix sub( A )=A(IA:IA+N-1,JA:JA+N-1).
  The factorization has the form
            sub( A ) = U' * U ,  if UPLO = 'U', or
            sub( A ) = L  * L',  if UPLO = 'L',
  where U is an upper triangular matrix and L is lower triangular.
  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 p x q.
  LOCr( K ) denotes the number of elements of K that a process
  would receive if K were distributed over the p 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 q 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
  This routine requires N <= NB_A-MOD(JA-1, NB_A) and square block
  decomposition ( MB_A = NB_A ).
  Arguments
  =========
  UPLO    (global input) CHARACTER
          = 'U':  Upper triangle of sub( A ) is stored;
          = 'L':  Lower triangle of sub( A ) is stored.
  N       (global input) INTEGER
          The number of rows and columns to be operated on, i.e. the
          order of the distributed submatrix sub( A ). N >= 0.
  A       (local input/local output) DOUBLE PRECISION pointer into the
          local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
          On entry, this array contains the local pieces of the
          N-by-N symmetric distributed matrix sub( A ) to be factored.
          If UPLO = 'U', the leading N-by-N upper triangular part of
          sub( A ) contains the upper triangular part of the matrix,
          and its strictly lower triangular part is not referenced.
          If UPLO = 'L', the leading N-by-N lower triangular part of
          sub( A ) contains the lower triangular part of the distribu-
          ted matrix, and its strictly upper triangular part is not
          referenced.  On exit, if UPLO = 'U', the upper triangular
          part of the distributed matrix contains the Cholesky factor
          U, if UPLO = 'L', the lower triangular part of the distribu-
          ted matrix contains the Cholesky factor L.
  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 and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix A.
  INFO    (local 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 = K, the leading minor of order K,
                A(IA:IA+K-1,JA:JA+K-1) is not positive definite, and
                the factorization could not be completed.
  =====================================================================
     .. Parameters ..

Display dynamic version

FontSize: - +

001        SUBROUTINE PDPOTF2( UPLO , N , A , IA , JA , DESCA , INFO )
002  
003  *     -- ScaLAPACK routine(version 1.7) --
004  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
005  *     and University of California , Berkeley.
006  *     May 1 , 1997
007  
008  *     .. Scalar Arguments ..
009        CHARACTER UPLO
010        INTEGER IA , INFO , JA , N
011        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
012       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
013        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
014       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
015       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
016        DOUBLE PRECISION ONE , ZERO
017        PARAMETER( ONE = 1.0D + 0 , ZERO = 0.0D + 0 )
018  *     ..
019  *     .. Local Scalars ..
020        LOGICAL UPPER
021        CHARACTER COLBTOP , ROWBTOP
022        INTEGER IACOL , IAROW , ICOFF , ICTXT , ICURR , IDIAG , IIA ,
023       $IOFFA , IROFF , J , JJA , LDA , MYCOL , MYROW ,
024       $NPCOL , NPROW
025        DOUBLE PRECISION AJJ
026  *     ..
027  *     .. External Subroutines ..
028        EXTERNAL BLACS_ABORT , BLACS_GRIDINFO , CHK1MAT , DGEMV ,
029       $DSCAL , IGEBR2D , IGEBS2D , INFOG2L , PB_TOPGET ,
030       $PXERBLA
031  *     ..
032  *     .. Intrinsic Functions ..
033        INTRINSIC MOD , SQRT
034  *     ..
035  *     .. External Functions ..
036        LOGICAL LSAME
037        DOUBLE PRECISION DDOT
038        EXTERNAL LSAME , DDOT
039  *     ..
040  *     .. Executable Statements ..
041  
042  *     Get grid parameters
043  
044        ICTXT = DESCA( CTXT_ )
045        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
046  
047  *     Test the input parameters.
048  
049        INFO = 0
050        IF( NPROW.EQ. - 1 ) THEN
050               
051            INFO = - (600 + CTXT_)
052        ELSE
052               
053            CALL CHK1MAT( N , 2 , N , 2 , IA , JA , DESCA , 6 , INFO )
054            IF( INFO.EQ.0 ) THEN
054                   
055                UPPER = LSAME( UPLO , 'U' )
056                IROFF = MOD( IA - 1 , DESCA( MB_ ) )
057                ICOFF = MOD( JA - 1 , DESCA( NB_ ) )
058                IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO , 'L' ) ) THEN
058                       
059                    INFO = - 1
060                ELSE IF( N + ICOFF.GT.DESCA( NB_ ) ) THEN
060                       
061                    INFO = - 2
062                ELSE IF( IROFF.NE.0 ) THEN
062                       
063                    INFO = - 4
064                ELSE IF( ICOFF.NE.0 ) THEN
064                       
065                    INFO = - 5
066                ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
066                       
067                    INFO = - (600 + NB_)
068                END IF
069            END IF
070        END IF
071  
072        IF( INFO.NE.0 ) THEN
072               
073            CALL PXERBLA( ICTXT , 'PDPOTF2' , - INFO )
074            CALL BLACS_ABORT( ICTXT , 1 )
075            RETURN
076        END IF
077  
078  *     Quick return if possible
079  
080        IF( N.EQ.0 )
080               
081       $    RETURN
082  
083  *         Compute local information
084  
085            CALL INFOG2L( IA , JA , DESCA , NPROW , NPCOL , MYROW , MYCOL , IIA , JJA ,
086       $    IAROW , IACOL )
087            CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Rowwise' , ROWBTOP )
088            CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Columnwise' , COLBTOP )
089  
090            IF( UPPER ) THEN
091  
092  *             Process(IAROW , IACOL) owns block to be factorized
093  
093                   
094                IF( MYROW.EQ.IAROW ) THEN
094                       
095                    IF( MYCOL.EQ.IACOL ) THEN
096  
097  *                     Compute the Cholesky factorization A = U'*U.
098  
098                           
099                        LDA = DESCA( LLD_ )
100                        IDIAG = IIA + ( JJA - 1 ) * LDA
101                        IOFFA = IDIAG
102  
103                        DO 10 J = JA , JA + N - 1
104  
105  *                         Compute U(J , J) and test for non - positive - definiteness.
106  
106                               
107                            AJJ = A( IDIAG ) -
108       $                    DDOT( J - JA , A( IOFFA ) , 1 , A( IOFFA ) , 1 )
109                            IF( AJJ.LE.ZERO ) THEN
109                                   
110                                A( IDIAG ) = AJJ
111                                INFO = J - JA + 1
112                                GO TO 20
113                            END IF
114                            AJJ = SQRT( AJJ )
115                            A( IDIAG ) = AJJ
116  
117  *                         Compute elements J + 1 : JA + N - 1 of row J.
118  
119                            IF( J.LT.JA + N - 1 ) THEN
119                                   
120                                ICURR = IDIAG + LDA
121                                CALL DGEMV( 'Transpose' , J - JA , JA + N - J - 1 , - ONE ,
122       $                        A( IOFFA + LDA ) , LDA , A( IOFFA ) , 1 ,
123       $                        ONE , A( ICURR ) , LDA )
124                                CALL DSCAL( N - J + JA - 1 , ONE / AJJ , A( ICURR ) , LDA )
125                            END IF
126                            IDIAG = IDIAG + LDA + 1
127                            IOFFA = IOFFA + LDA
128     10                 CONTINUE
129  
130     20 CONTINUE
131  
132  *     Broadcast INFO to all processes in my IAROW.
133  
134        CALL IGEBS2D( ICTXT , 'Rowwise' , ROWBTOP , 1 , 1 , INFO , 1 )
135  
136        ELSE
137  
137               
138            CALL IGEBR2D( ICTXT , 'Rowwise' , ROWBTOP , 1 , 1 , INFO , 1 ,
139       $    MYROW , IACOL )
140        END IF
141  
142  *     IAROW bcasts along columns so that everyone has INFO
143  
144        CALL IGEBS2D( ICTXT , 'Columnwise' , COLBTOP , 1 , 1 , INFO , 1 )
145  
146        ELSE
147  
147               
148            CALL IGEBR2D( ICTXT , 'Columnwise' , COLBTOP , 1 , 1 , INFO , 1 ,
149       $    IAROW , MYCOL )
150  
151        END IF
152  
153        ELSE
154  
155  *         Process(IAROW , IACOL) owns block to be factorized
156  
156               
157            IF( MYCOL.EQ.IACOL ) THEN
157                   
158                IF( MYROW.EQ.IAROW ) THEN
159  
160  *                 Compute the Cholesky factorization A = L*L'.
161  
161                       
162                    LDA = DESCA( LLD_ )
163                    IDIAG = IIA + ( JJA - 1 ) * LDA
164                    IOFFA = IDIAG
165  
166                    DO 30 J = JA , JA + N - 1
167  
168  *                     Compute L(J , J) and test for non - positive - definiteness.
169  
169                           
170                        AJJ = A( IDIAG ) -
171       $                DDOT( J - JA , A( IOFFA ) , LDA , A( IOFFA ) , LDA )
172                        IF( AJJ.LE.ZERO ) THEN
172                               
173                            A( IDIAG ) = AJJ
174                            INFO = J - JA + 1
175                            GO TO 40
176                        END IF
177                        AJJ = SQRT( AJJ )
178                        A( IDIAG ) = AJJ
179  
180  *                     Compute elements J + 1 : JA + N - 1 of column J.
181  
182                        IF( J.LT.JA + N - 1 ) THEN
182                               
183                            ICURR = IDIAG + 1
184                            CALL DGEMV( 'No transpose' , JA + N - J - 1 , J - JA , - ONE ,
185       $                    A( IOFFA + 1 ) , LDA , A( IOFFA ) , LDA ,
186       $                    ONE , A( ICURR ) , 1 )
187                            CALL DSCAL( JA + N - J - 1 , ONE / AJJ , A( ICURR ) , 1 )
188                        END IF
189                        IDIAG = IDIAG + LDA + 1
190                        IOFFA = IOFFA + 1
191     30             CONTINUE
192  
193     40 CONTINUE
194  
195  *     Broadcast INFO to everyone in IACOL
196  
197        CALL IGEBS2D( ICTXT , 'Columnwise' , COLBTOP , 1 , 1 , INFO ,
198       $1 )
199  
200        ELSE
201  
201               
202            CALL IGEBR2D( ICTXT , 'Columnwise' , COLBTOP , 1 , 1 , INFO ,
203       $    1 , IAROW , MYCOL )
204  
205        END IF
206  
207  *     IACOL bcasts INFO along rows so that everyone has it
208  
209        CALL IGEBS2D( ICTXT , 'Rowwise' , ROWBTOP , 1 , 1 , INFO , 1 )
210  
211        ELSE
212  
212               
213            CALL IGEBR2D( ICTXT , 'Rowwise' , ROWBTOP , 1 , 1 , INFO , 1 ,
214       $    MYROW , IACOL )
215  
216        END IF
217  
218        END IF
219  
220        RETURN
221  
222  *     End of PDPOTF2
223  
224        END3256

Variables in Routine PDPOTF2()

Summary Report
Data Type	Quantity	Size(byte)
CHARACTER	3	3
DOUBLE PRECISION	4	16
INTEGER	31	124
LOGICAL	2	2
REAL	1	4
TOTAL	41	149

List of Variables

CHARACTER

COLBTOP ROWBTOP UPLO

DOUBLE PRECISION

AJJ DDOT ONE ZERO

INTEGER

BLOCK_CYCLIC_2D CSRC_ CTXT_ DLEN_ DTYPE_
IA IACOL IAROW ICOFF ICTXT
ICURR IDIAG IIA INFO IOFFA
IROFF J JA JJA LDA
LLD_ M_ MB_ MYCOL MYROW
N N_ NB_ NPCOL NPROW
RSRC_

LOGICAL

LSAME UPPER

REAL

A

Variables Dependence Graph
Put the mouse over a right hand side variable to display the corresponding line of the dependence
-
- -
A <--- AJJA( IDIAG ) = AJJ{2A( IDIAG ) = AJJ, 3A( IDIAG ) = AJJ, 4A( IDIAG ) = AJJ}

AJJ <--- AAJJ = A( IDIAG ) -{2AJJ = A( IDIAG ) -}, IDIAGAJJ = A( IDIAG ) -{2AJJ = A( IDIAG ) -}, IOFFAAJJ = A( IDIAG ) -{2AJJ = A( IDIAG ) -}, AJJAJJ = SQRT( AJJ ){2AJJ = SQRT( AJJ )}, JAJJ = A( IDIAG ) -{2AJJ = A( IDIAG ) -}, JAAJJ = A( IDIAG ) -{2AJJ = A( IDIAG ) -}, LDAAJJ = A( IDIAG ) -, DDOTAJJ = A( IDIAG ) -{2AJJ = A( IDIAG ) -}

ICOFF <--- JAICOFF = MOD( JA-1, DESCA( NB_ ) ), NB_ICOFF = MOD( JA-1, DESCA( NB_ ) )

ICTXT <--- CTXT_ICTXT = DESCA( CTXT_ )

ICURR <--- IDIAGICURR = IDIAG + LDA{2ICURR = IDIAG + 1}, LDAICURR = IDIAG + LDA

IDIAG <--- IDIAGIDIAG = IDIAG + LDA + 1{2IDIAG = IDIAG + LDA + 1}, IIAIDIAG = IIA + ( JJA - 1 ) * LDA{2IDIAG = IIA + ( JJA - 1 ) * LDA}, JJAIDIAG = IIA + ( JJA - 1 ) * LDA{2IDIAG = IIA + ( JJA - 1 ) * LDA}, LDAIDIAG = IIA + ( JJA - 1 ) * LDA{2IDIAG = IDIAG + LDA + 1, 3IDIAG = IIA + ( JJA - 1 ) * LDA, 4IDIAG = IDIAG + LDA + 1}

INFO <--- JINFO = J - JA + 1{2INFO = J - JA + 1}, JAINFO = J - JA + 1{2INFO = J - JA + 1}, NB_INFO = -(600+NB_), CTXT_INFO = -(600+CTXT_)

IOFFA <--- IDIAGIOFFA = IDIAG{2IOFFA = IDIAG}, IOFFAIOFFA = IOFFA + LDA{2IOFFA = IOFFA + 1}, LDAIOFFA = IOFFA + LDA

IROFF <--- IAIROFF = MOD( IA-1, DESCA( MB_ ) ), MB_IROFF = MOD( IA-1, DESCA( MB_ ) )

J <--- JADO 10 J = JA, JA+N-1{2DO 30 J = JA, JA+N-1}, NDO 10 J = JA, JA+N-1{2DO 30 J = JA, JA+N-1}

LDA <--- LLD_LDA = DESCA( LLD_ ){2LDA = DESCA( LLD_ )}

UPPER <--- LSAMEUPPER = LSAME( UPLO, 'U' ), UPLOUPPER = LSAME( UPLO, 'U' )