|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 REAL FUNCTION PSQRT14( TRANS, M, N, NRHS, A, IA, JA, 00002 $ DESCA, X, IX, JX, DESCX, WORK ) 00003 * 00004 * -- ScaLAPACK routine (version 1.7) -- 00005 * University of Tennessee, Knoxville, Oak Ridge National Laboratory, 00006 * and University of California, Berkeley. 00007 * May 1, 1997 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER TRANS 00011 INTEGER IA, IX, JA, JX, M, N, NRHS 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER DESCA( * ), DESCX( * ) 00015 REAL A( * ), WORK( * ), X( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * PSQRT14 checks whether sub( X ) is in the row space of sub( A ) or 00022 * sub( A )', where sub( A ) denotes A( IA:IA+M-1, JA:JA+N-1 ) and 00023 * sub( X ) denotes X( IX:IX+N-1, JX:JX+NRHS-1 ) if TRANS = 'N', and 00024 * X( IX:IX+N-1, JX:JX+NRHS-1 ) otherwise. It does so by scaling both 00025 * sub( X ) and sub( A ) such that their norms are in the range 00026 * [sqrt(eps), 1/sqrt(eps)], then computing an LQ factorization of 00027 * [sub( A )',sub( X )]' (if TRANS = 'N') or a QR factorization of 00028 * [sub( A ),sub( X )] otherwise, and returning the norm of the trailing 00029 * triangle, scaled by MAX(M,N,NRHS)*eps. 00030 * 00031 * Notes 00032 * ===== 00033 * 00034 * Each global data object is described by an associated description 00035 * vector. This vector stores the information required to establish 00036 * the mapping between an object element and its corresponding process 00037 * and memory location. 00038 * 00039 * Let A be a generic term for any 2D block cyclicly distributed array. 00040 * Such a global array has an associated description vector DESCA. 00041 * In the following comments, the character _ should be read as 00042 * "of the global array". 00043 * 00044 * NOTATION STORED IN EXPLANATION 00045 * --------------- -------------- -------------------------------------- 00046 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case, 00047 * DTYPE_A = 1. 00048 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating 00049 * the BLACS process grid A is distribu- 00050 * ted over. The context itself is glo- 00051 * bal, but the handle (the integer 00052 * value) may vary. 00053 * M_A (global) DESCA( M_ ) The number of rows in the global 00054 * array A. 00055 * N_A (global) DESCA( N_ ) The number of columns in the global 00056 * array A. 00057 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute 00058 * the rows of the array. 00059 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute 00060 * the columns of the array. 00061 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first 00062 * row of the array A is distributed. 00063 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the 00064 * first column of the array A is 00065 * distributed. 00066 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local 00067 * array. LLD_A >= MAX(1,LOCr(M_A)). 00068 * 00069 * Let K be the number of rows or columns of a distributed matrix, 00070 * and assume that its process grid has dimension p x q. 00071 * LOCr( K ) denotes the number of elements of K that a process 00072 * would receive if K were distributed over the p processes of its 00073 * process column. 00074 * Similarly, LOCc( K ) denotes the number of elements of K that a 00075 * process would receive if K were distributed over the q processes of 00076 * its process row. 00077 * The values of LOCr() and LOCc() may be determined via a call to the 00078 * ScaLAPACK tool function, NUMROC: 00079 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ), 00080 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ). 00081 * An upper bound for these quantities may be computed by: 00082 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A 00083 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A 00084 * 00085 * Arguments 00086 * ========= 00087 * 00088 * TRANS (global input) CHARACTER*1 00089 * = 'N': No transpose, check for sub( X ) in the row space of 00090 * sub( A ), 00091 * = 'T': Transpose, check for sub( X ) in row space of 00092 * sub( A )'. 00093 * 00094 * M (global input) INTEGER 00095 * The number of rows to be operated on, i.e. the number of rows 00096 * of the distributed submatrix sub( A ). M >= 0. 00097 * 00098 * N (global input) INTEGER 00099 * The number of columns to be operated on, i.e. the number of 00100 * columns of the distributed submatrix sub( A ). N >= 0. 00101 * 00102 * NRHS (global input) INTEGER 00103 * The number of right hand sides, i.e., the number of columns 00104 * of the distributed submatrix sub( X ). NRHS >= 0. 00105 * 00106 * A (local input) REAL pointer into the local memory 00107 * to an array of dimension (LLD_A, LOCc(JA+N-1)). This array 00108 * contains the local pieces of the M-by-N distributed matrix 00109 * sub( A ). 00110 * 00111 * IA (global input) INTEGER 00112 * The row index in the global array A indicating the first 00113 * row of sub( A ). 00114 * 00115 * JA (global input) INTEGER 00116 * The column index in the global array A indicating the 00117 * first column of sub( A ). 00118 * 00119 * DESCA (global and local input) INTEGER array of dimension DLEN_. 00120 * The array descriptor for the distributed matrix A. 00121 * 00122 * X (local input) REAL pointer into the local 00123 * memory to an array of dimension (LLD_X,LOCc(JX+NRHS-1)). 00124 * On entry, this array contains the local pieces of the 00125 * N-by-NRHS distributed submatrix sub( X ) if TRANS = 'N', 00126 * and the M-by-NRHS distributed submatrix sub( X ) otherwise. 00127 * 00128 * IX (global input) INTEGER 00129 * The row index in the global array X indicating the first 00130 * row of sub( X ). 00131 * 00132 * JX (global input) INTEGER 00133 * The column index in the global array X indicating the 00134 * first column of sub( X ). 00135 * 00136 * DESCX (global and local input) INTEGER array of dimension DLEN_. 00137 * The array descriptor for the distributed matrix X. 00138 * 00139 * WORK (local workspace) REAL array dimension (LWORK) 00140 * If TRANS='N', LWORK >= MNRHSP * NQ + LTAU + LWF and 00141 * LWORK >= MP * NNRHSQ + LTAU + LWF otherwise, where 00142 * 00143 * IF TRANS='N', (LQ fact) 00144 * MNRHSP = NUMROC( M+NRHS+IROFFA, MB_A, MYROW, IAROW, 00145 * NPROW ) 00146 * LTAU = NUMROC( IA+MIN( M+NRHS, N )-1, MB_A, MYROW, 00147 * RSRC_A, NPROW ) 00148 * LWF = MB_A * ( MB_A + MNRHSP + NQ0 ) 00149 * ELSE (QR fact) 00150 * NNRHSQ = NUMROC( N+NRHS+ICOFFA, NB_A, MYCOL, IACOL, 00151 * NPCOL ) 00152 * LTAU = NUMROC( JA+MIN( M, N+NRHS )-1, NB_A, MYCOL, 00153 * CSRC_A, NPCOL ) 00154 * LWF = NB_A * ( NB_A + MP0 + NNRHSQ ) 00155 * END IF 00156 * 00157 * and, 00158 * 00159 * IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ), 00160 * IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ), 00161 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ), 00162 * MP0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ), 00163 * NQ0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ). 00164 * 00165 * INDXG2P and NUMROC are ScaLAPACK tool functions; 00166 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling 00167 * the subroutine BLACS_GRIDINFO. 00168 * 00169 * 00170 * ===================================================================== 00171 * 00172 * .. Parameters .. 00173 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 00174 $ LLD_, MB_, M_, NB_, N_, RSRC_ 00175 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00176 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00177 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00178 REAL ONE, ZERO 00179 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00180 * .. 00181 * .. Local Scalars .. 00182 LOGICAL TPSD 00183 INTEGER IACOL, IAROW, ICOFFA, ICTXT, IDUM, IIA, INFO, 00184 $ IPTAU, IPW, IPWA, IROFFA, IWA, IWX, J, JJA, 00185 $ JWA, JWX, LDW, LWORK, MPWA, MPW, MQW, MYCOL, 00186 $ MYROW, NPCOL, NPROW, NPW, NQWA, NQW 00187 REAL AMAX, ANRM, ERR, XNRM 00188 * .. 00189 * .. Local Arrays .. 00190 INTEGER DESCW( DLEN_ ), IDUM1( 1 ), IDUM2( 1 ) 00191 REAL RWORK( 1 ) 00192 * .. 00193 * .. External Functions .. 00194 LOGICAL LSAME 00195 INTEGER NUMROC 00196 REAL PSLANGE, PSLAMCH 00197 EXTERNAL LSAME, NUMROC, PSLANGE, PSLAMCH 00198 * .. 00199 * .. External Subroutines .. 00200 EXTERNAL BLACS_GRIDINFO, DESCSET, INFOG2L, PSAMAX, 00201 $ PSCOPY, PSGELQF, PSGEQRF, PSLACPY, 00202 $ PSLASCL, PXERBLA, SGAMX2D 00203 * .. 00204 * .. Intrinsic Functions .. 00205 INTRINSIC ABS, MAX, MIN, MOD, REAL 00206 * .. 00207 * .. Executable Statements .. 00208 * 00209 * Get grid parameters 00210 * 00211 ICTXT = DESCA( CTXT_ ) 00212 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00213 * 00214 PSQRT14 = ZERO 00215 * 00216 IPWA = 1 00217 IROFFA = MOD( IA-1, DESCA( MB_ ) ) 00218 ICOFFA = MOD( JA-1, DESCA( NB_ ) ) 00219 IWA = IROFFA + 1 00220 JWA = ICOFFA + 1 00221 CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, 00222 $ JJA, IAROW, IACOL ) 00223 MPWA = NUMROC( M+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW ) 00224 NQWA = NUMROC( N+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 00225 * 00226 INFO = 0 00227 IF( LSAME( TRANS, 'N' ) ) THEN 00228 IF( N.LE.0 .OR. NRHS.LE.0 ) 00229 $ RETURN 00230 TPSD = .FALSE. 00231 MPW = NUMROC( M+NRHS+IROFFA, DESCA( MB_ ), MYROW, IAROW, 00232 $ NPROW ) 00233 NQW = NQWA 00234 * 00235 * Assign descriptor DESCW for workspace WORK and pointers to 00236 * matrices sub( A ) and sub( X ) in workspace 00237 * 00238 IWX = IWA + M 00239 JWX = JWA 00240 LDW = MAX( 1, MPW ) 00241 CALL DESCSET( DESCW, M+NRHS+IROFFA, N+ICOFFA, DESCA( MB_ ), 00242 $ DESCA( NB_ ), IAROW, IACOL, ICTXT, LDW ) 00243 * 00244 ELSE IF( LSAME( TRANS, 'T' ) ) THEN 00245 IF( M.LE.0 .OR. NRHS.LE.0 ) 00246 $ RETURN 00247 TPSD = .TRUE. 00248 MPW = MPWA 00249 NQW = NUMROC( N+NRHS+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, 00250 $ NPCOL ) 00251 * 00252 * Assign descriptor DESCW for workspace WORK and pointers to 00253 * matrices sub( A ) and sub( X ) in workspace 00254 * 00255 IWX = IWA 00256 JWX = JWA + N 00257 LDW = MAX( 1, MPW ) 00258 CALL DESCSET( DESCW, M+IROFFA, N+NRHS+ICOFFA, DESCA( MB_ ), 00259 $ DESCA( NB_ ), IAROW, IACOL, ICTXT, LDW ) 00260 ELSE 00261 CALL PXERBLA( ICTXT, 'PSQRT14', -1 ) 00262 RETURN 00263 END IF 00264 * 00265 * Copy and scale sub( A ) 00266 * 00267 IPTAU = IPWA + MPW*NQW 00268 CALL PSLACPY( 'All', M, N, A, IA, JA, DESCA, WORK( IPWA ), IWA, 00269 $ JWA, DESCW ) 00270 RWORK( 1 ) = ZERO 00271 ANRM = PSLANGE( 'M', M, N, WORK( IPWA ), IWA, JWA, DESCW, RWORK ) 00272 IF( ANRM.NE.ZERO ) 00273 $ CALL PSLASCL( 'G', ANRM, ONE, M, N, WORK( IPWA ), IWA, 00274 $ JWA, DESCW, INFO ) 00275 * 00276 * Copy sub( X ) or sub( X )' into the right place and scale it 00277 * 00278 IF( TPSD ) THEN 00279 * 00280 * Copy sub( X ) into columns jwa+n:jwa+n+nrhs-1 of work 00281 * 00282 DO 10 J = 1, NRHS 00283 CALL PSCOPY( M, X, IX, JX+J-1, DESCX, 1, WORK( IPWA ), IWX, 00284 $ JWX+J-1, DESCW, 1 ) 00285 10 CONTINUE 00286 XNRM = PSLANGE( 'M', M, NRHS, WORK( IPWA ), IWX, JWX, DESCW, 00287 $ RWORK ) 00288 IF( XNRM.NE.ZERO ) 00289 $ CALL PSLASCL( 'G', XNRM, ONE, M, NRHS, WORK( IPWA ), IWX, 00290 $ JWX, DESCW, INFO ) 00291 * 00292 * Compute QR factorization of work(iwa:iwa+m-1,jwa:jwa+n+nrhs-1) 00293 * 00294 MQW = NUMROC( M+ICOFFA, DESCA( NB_ ), MYCOL, IACOL, NPCOL ) 00295 IPW = IPTAU + MIN( MQW, NQW ) 00296 LWORK = DESCW( NB_ ) * ( MPW + NQW + DESCW( NB_ ) ) 00297 CALL PSGEQRF( M, N+NRHS, WORK( IPWA ), IWA, JWA, DESCW, 00298 $ WORK( IPTAU ), WORK( IPW ), LWORK, INFO ) 00299 * 00300 * Compute largest entry in upper triangle of 00301 * work(iwa+n:iwa+m-1,jwa+n:jwa+n+nrhs-1) 00302 * 00303 ERR = ZERO 00304 IF( N.LT.M ) THEN 00305 DO 20 J = JWX, JWA+N+NRHS-1 00306 CALL PSAMAX( MIN(M-N,J-JWX+1), AMAX, IDUM, WORK( IPWA ), 00307 $ IWA+N, J, DESCW, 1 ) 00308 ERR = MAX( ERR, ABS( AMAX ) ) 00309 20 CONTINUE 00310 END IF 00311 CALL SGAMX2D( ICTXT, 'All', ' ', 1, 1, ERR, 1, IDUM1, IDUM2, 00312 $ -1, -1, 0 ) 00313 * 00314 ELSE 00315 * 00316 * Copy sub( X )' into rows iwa+m:iwa+m+nrhs-1 of work 00317 * 00318 DO 30 J = 1, NRHS 00319 CALL PSCOPY( N, X, IX, JX+J-1, DESCX, 1, WORK( IPWA ), 00320 $ IWX+J-1, JWX, DESCW, DESCW( M_ ) ) 00321 30 CONTINUE 00322 * 00323 XNRM = PSLANGE( 'M', NRHS, N, WORK( IPWA ), IWX, JWX, DESCW, 00324 $ RWORK ) 00325 IF( XNRM.NE.ZERO ) 00326 $ CALL PSLASCL( 'G', XNRM, ONE, NRHS, N, WORK( IPWA ), IWX, 00327 $ JWX, DESCW, INFO ) 00328 * 00329 * Compute LQ factorization of work(iwa:iwa+m+nrhs-1,jwa:jwa+n-1) 00330 * 00331 NPW = NUMROC( N+IROFFA, DESCA( MB_ ), MYROW, IAROW, NPROW ) 00332 IPW = IPTAU + MIN( MPW, NPW ) 00333 LWORK = DESCW( MB_ ) * ( MPW + NQW + DESCW( MB_ ) ) 00334 CALL PSGELQF( M+NRHS, N, WORK( IPWA ), IWA, JWA, DESCW, 00335 $ WORK( IPTAU ), WORK( IPW ), LWORK, INFO ) 00336 * 00337 * Compute largest entry in lower triangle in 00338 * work(iwa+m:iwa+m+nrhs-1,jwa+m:jwa+n-1) 00339 * 00340 ERR = ZERO 00341 DO 40 J = JWA+M, MIN( JWA+N-1, JWA+M+NRHS-1 ) 00342 CALL PSAMAX( JWA+M+NRHS-J, AMAX, IDUM, WORK( IPWA ), 00343 $ IWX+J-JWA-M, J, DESCW, 1 ) 00344 ERR = MAX( ERR, ABS( AMAX ) ) 00345 40 CONTINUE 00346 CALL SGAMX2D( ICTXT, 'All', ' ', 1, 1, ERR, 1, IDUM1, IDUM2, 00347 $ -1, -1, 0 ) 00348 * 00349 END IF 00350 * 00351 PSQRT14 = ERR / ( REAL( MAX( M, N, NRHS ) ) * 00352 $ PSLAMCH( ICTXT, 'Epsilon' ) ) 00353 * 00354 RETURN 00355 * 00356 * End of PSQRT14 00357 * 00358 END
