Routine: PCPOTF2()  File: SRC\pcpotf2.f

 
 
# lines: 358
  # code: 358
  # comment: 0
  # blank:0
# Variables:42
# Callers:1
# Callings:0
# Words:190
# Keywords:122
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PCPOTF2 computes the Cholesky factorization of a complex hermitian
  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) COMPLEX 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 Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

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