ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pdsvdchk.f
Go to the documentation of this file.
00001       SUBROUTINE PDSVDCHK( M, N, A, IA, JA, DESCA, U, IU, JU, DESCU, VT,
00002      $                     IVT, JVT, DESCVT, S, THRESH, WORK, LWORK,
00003      $                     RESULT, CHK, MTM )
00004 *
00005 *  -- ScaLAPACK routine (version 1.7) --
00006 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00007 *     and University of California, Berkeley.
00008 *     May 1, 1997
00009 *
00010 *     .. Scalar Arguments ..
00011       INTEGER            IA, IU, IVT, JA, JU, JVT, LWORK, M, N
00012       DOUBLE PRECISION   CHK, MTM, THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       INTEGER            DESCA( * ), DESCU( * ), DESCVT( * ),
00016      $                   RESULT( * )
00017       DOUBLE PRECISION   A( * ), S( * ), U( * ), VT( * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  For given two-dimensional matrices A, U, VT, and one-dimensional
00024 *  array D compute the following four tests:
00025 *
00026 *  (1)   | A - U*diag(S) VT | / ( |A| max(M,N) ulp )
00027 *
00028 *  (2)   | I - U'*U | / ( M ulp )
00029 *
00030 *  (3)   | I - VT*VT' | / ( N ulp ),
00031 *
00032 *  (4)   S contains SIZE = MIN( M, N )  nonnegative values in
00033 *   decreasing order.
00034 *   It then compares result of computations (1)-(3)
00035 *   with TRESH and returns results of comparisons and test (4) in
00036 *   RESULT(I). When the i-th test fails, value of RESULT( I ) is set
00037 *   to 1.
00038 *
00039 *  Notes
00040 *  =====
00041 *
00042 *  Each global data object is described by an associated description
00043 *  vector.  This vector stores the information required to establish
00044 *  the mapping between an object element and its corresponding process
00045 *  and memory location.
00046 *
00047 *  Let A be a generic term for any 2D block cyclicly distributed array.
00048 *  Such a global array has an associated description vector DESCA.
00049 *  In the following comments, the character _ should be read as
00050 *  "of the global array".
00051 *
00052 *  NOTATION        STORED IN      EXPLANATION
00053 *  --------------- -------------- --------------------------------------
00054 *  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
00055 *                                 DTYPE_A = 1.
00056 *  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
00057 *                                 the BLACS process grid A is distribu-
00058 *                                 ted over. The context itself is glo-
00059 *                                 bal, but the handle (the integer
00060 *                                 value) may vary.
00061 *  M_A    (global) DESCA( M_ )    The number of rows in the global
00062 *                                 array A.
00063 *  N_A    (global) DESCA( N_ )    The number of columns in the global
00064 *                                 array A.
00065 *  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
00066 *                                 the rows of the array.
00067 *  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
00068 *                                 the columns of the array.
00069 *  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
00070 *                                 row of the array A is distributed.
00071 *  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
00072 *                                 first column of the array A is
00073 *                                 distributed.
00074 *  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
00075 *                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
00076 *
00077 *  Let K be the number of rows or columns of a distributed matrix,
00078 *  and assume that its process grid has dimension p x q.
00079 *  LOCr( K ) denotes the number of elements of K that a process
00080 *  would receive if K were distributed over the p processes of its
00081 *  process column.
00082 *  Similarly, LOCc( K ) denotes the number of elements of K that a
00083 *  process would receive if K were distributed over the q processes of
00084 *  its process row.
00085 *  The values of LOCr() and LOCc() may be determined via a call to the
00086 *  ScaLAPACK tool function, NUMROC:
00087 *          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
00088 *          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
00089 *  An upper bound for these quantities may be computed by:
00090 *          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
00091 *          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
00092 *
00093 *  Arguments
00094 *  =========
00095 *
00096 *     MP = number of local rows in A and  U
00097 *     NQ = number of local columns in A and VT
00098 *     SIZEP = number of local rows in VT
00099 *     SIZEQ = number of local columns in U
00100 *
00101 *  M      (global input) INTEGER
00102 *          Matrix size.
00103 *          The number of global rows in A and U and
00104 *
00105 *  N      (global input) INTEGER
00106 *          The number of global columns in A and VT.
00107 *
00108 *  A       (input) block cyclic distributed DOUBLE PRECISION array,
00109 *          global dimension (M, N), local dimension (DESCA( DLEN_ ), NQ)
00110 *          Contains the original test matrix.
00111 *
00112 *  IA      (global input) INTEGER
00113 *          The global row index of the submatrix of the distributed
00114 *          matrix A to operate on.
00115 *
00116 *  JA      (global input) INTEGER
00117 *          The global column index of the submatrix of the distributed
00118 *          matrix A to operate on.
00119 *
00120 *  DESCA   (global and local input) INTEGER array of dimension DLEN_
00121 *          The array descriptor for the distributed matrix A.
00122 *
00123 *  U       (local input) DOUBLE PRECISION array
00124 *           global dimension (M, SIZE), local dimension
00125 *          (DESCU( DLEN_ ), SIZEQ)
00126 *           Contains left singular vectors of matrix A.
00127 *
00128 *  IU      (global input) INTEGER
00129 *          The global row index of the submatrix of the distributed
00130 *          matrix U to operate on.
00131 *
00132 *  JU      (global input) INTEGER
00133 *          The global column index of the submatrix of the distributed
00134 *          matrix U to operate on.
00135 *
00136 *  DESCU   (global and local input) INTEGER array of dimension DLEN_
00137 *          The array descriptor for the distributed matrix U.
00138 *
00139 *  VT       (local input) DOUBLE PRECISION array
00140 *           global dimension (SIZE, N), local dimension
00141 *           (DESCVT( DLEN_ ), NQ)
00142 *           Contains right singular vectors of matrix A.
00143 *
00144 *  IVT     (global input) INTEGER
00145 *          The global row index of the submatrix of the distributed
00146 *          matrix VT to operate on.
00147 *
00148 *  JVT      (global input) INTEGER
00149 *          The global column index of the submatrix of the distributed
00150 *          matrix VT to operate on.
00151 *
00152 *  DESCVT   (global and local input) INTEGER array of dimension DLEN_
00153 *          The array descriptor for the distributed matrix VT.
00154 *
00155 *  S       (global input) DOUBLE PRECISION array, dimension (SIZE)
00156 *          Contains the computed singular values
00157 *
00158 *  THRESH  (input) DOUBLE PRECISION
00159 *          A test will count as "failed" if the "error", computed as
00160 *          described below, exceeds THRESH.  Note that the error
00161 *          is scaled to be O(1), so THRESH should be a reasonably
00162 *          small multiple of 1, e.g., 10 or 100.  In particular,
00163 *          it should not depend on the precision (single vs. double)
00164 *          or the size of the matrix.  It must be at least zero.
00165 *
00166 *  WORK    (local workspace) DOUBLE PRECISION array, dimension (LWORK)
00167 *
00168 *  LWORK   (local input) INTEGER
00169 *          The length of the array WORK.
00170 *          LWORK >= 1 + SIZEQ*SIZEP + MAX[WORK(pdlange(size,size)),
00171 *          WORK(pdlange(m,n))],
00172 *          where
00173 *          SIZEQ = NUMROC( SIZE, DESCU( NB_ ), MYCOL, 0, NPCOL ),
00174 *          SIZEP = NUMROC( SIZE, DESCVT( MB_ ), MYROW, 0, NPROW ),
00175 *          and worekspaces required to call pdlange are
00176 *          WORK(pdlange(size,size)) < MAX(SIZEQ0,2) < SIZEB +2,
00177 *          WORK(pdlange(m,n)) < MAX(NQ0,2) < SIZEB +2,
00178 *          SIZEB = MAX(M, N)
00179 *          Finally, upper limit on required workspace is
00180 *          LWORK >  1 + SIZEQ*SIZEP + SIZEB + 2
00181 *
00182 *  RESULT  (global input/output) INTEGER array. Four first elements of
00183 *          the array are set to 0 or 1 depending on passing four
00184 *          respective tests ( see above in Purpose ). The elements of
00185 *          RESULT are set to
00186 *          0 if the test passes i.e.
00187 *            | A - U*diag(S)*VT | / ( |A| max(M,N) ulp ) <= THRESH
00188 *          1 if the test fails  i.e.
00189 *            | A - U*diag(S)*VT | / ( |A| max(M,N) ulp ) >  THRESH
00190 *
00191 *  CHK     (global output) DOUBLE PRECISION
00192 *           value of the | A - U*diag(S) VT | / ( |A| max(M,N) ulp )
00193 *
00194 *  MTM     (global output) DOUBLE PRECISION
00195 *           maximum of the two values:
00196 *            | I - U'*U | / ( M ulp ) and  | I - VT*VT' | / ( N ulp )
00197 *
00198 * ======================================================================
00199 *
00200 *     .. Parameters ..
00201       INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
00202      $                   MB_, NB_, RSRC_, CSRC_, LLD_
00203       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00204      $                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00205      $                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00206       DOUBLE PRECISION   ZERO, ONE, MONE
00207       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, MONE = -1.0D0 )
00208 *     ..
00209 *     .. Local Scalars ..
00210       INTEGER            I, INFO, LDR, LOCALCOL, LWMIN, MP, MX, MYCOL,
00211      $                   MYROW, NPCOL, NPROW, NQ, PCOL, PTRR, PTRWORK,
00212      $                   SIZE, SIZEP, SIZEPOS, SIZEQ
00213       DOUBLE PRECISION   FIRST, NORMA, NORMAI, NORMU, NORMVT, SECOND,
00214      $                   THRESHA, ULP
00215 *     ..
00216 *     .. Local Arrays ..
00217       INTEGER            DESCR( DLEN_ )
00218 *     ..
00219 *     .. External Functions ..
00220       INTEGER            INDXG2L, INDXG2P, NUMROC
00221       DOUBLE PRECISION   PDLAMCH, PDLANGE
00222       EXTERNAL           INDXG2L, INDXG2P, NUMROC, PDLAMCH, PDLANGE
00223 *     ..
00224 *     .. External Subroutines ..
00225       EXTERNAL           BLACS_GRIDINFO, CHK1MAT, DESCINIT, DSCAL,
00226      $                   PDELSET, PDGEMM, PDLASET, PXERBLA
00227 *     ..
00228 *     .. Intrinsic Functions ..
00229       INTRINSIC          MAX, MIN
00230 *     ..
00231 *     .. Executable Statements ..
00232 *     This is just to keep ftnchek happy
00233       IF( BLOCK_CYCLIC_2D*CSRC_*DTYPE_*M_*N_*RSRC_.LT.0 ) RETURN
00234 *
00235 *     Test the input parameters.
00236 *
00237       CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL )
00238       INFO = 0
00239       SIZE = MIN( M, N )
00240 *
00241 *     Sizepos is a number of parameters to pdsvdchk plus one. It's used
00242 *     for the error reporting.
00243 *
00244       SIZEPOS = 22
00245       IF( NPROW.EQ.-1 ) THEN
00246          INFO = -607
00247       ELSE
00248          CALL CHK1MAT( M, 1, N, 2, IA, JA, DESCA, 6, INFO )
00249          CALL CHK1MAT( M, 1, SIZE, SIZEPOS, IU, JU, DESCU, 10, INFO )
00250          CALL CHK1MAT( SIZE, SIZEPOS, N, 2, IVT, JVT, DESCVT, 14, INFO )
00251       END IF
00252 *
00253       IF( INFO.EQ.0 ) THEN
00254 *
00255 *     Calculate workspace
00256 *
00257          MP = NUMROC( M, DESCA( MB_ ), MYROW, 0, NPROW )
00258          NQ = NUMROC( N, DESCA( NB_ ), MYCOL, 0, NPCOL )
00259          SIZEP = NUMROC( SIZE, DESCVT( MB_ ), MYROW, 0, NPROW )
00260          SIZEQ = NUMROC( SIZE, DESCU( NB_ ), MYCOL, 0, NPCOL )
00261          MX = MAX( SIZEQ, NQ )
00262          LWMIN = 2 + SIZEQ*SIZEP + MAX( 2, MX )
00263          WORK( 1 ) = LWMIN
00264          IF( LWORK.EQ.-1 )
00265      $      GO TO 40
00266          IF( LWORK.LT.LWMIN ) THEN
00267             INFO = -18
00268          ELSE IF( THRESH.LE.0 ) THEN
00269             INFO = -16
00270          END IF
00271       END IF
00272       IF( INFO.NE.0 ) THEN
00273          CALL PXERBLA( DESCA( CTXT_ ), 'PDSVDCHK', -INFO )
00274          RETURN
00275       END IF
00276 *
00277       LDR = MAX( 1, SIZEP )
00278       ULP = PDLAMCH( DESCA( CTXT_ ), 'P' )
00279       NORMAI = PDLANGE( '1', M, N, A, IA, JA, DESCA, WORK )
00280 *
00281 *     Allocate array R of global dimension SIZE x SIZE for testing
00282 *
00283       PTRR = 2
00284       PTRWORK = PTRR + SIZEQ*SIZEP
00285 *
00286       CALL DESCINIT( DESCR, SIZE, SIZE, DESCVT( MB_ ), DESCU( NB_ ), 0,
00287      $               0, DESCA( CTXT_ ), LDR, INFO )
00288 *
00289 *     Test 2. Form identity matrix R  and make  check norm(U'*U - I )
00290 *
00291       CALL PDLASET( 'Full', SIZE, SIZE, ZERO, ONE, WORK( PTRR ), 1, 1,
00292      $              DESCR )
00293       CALL PDGEMM( 'T', 'N', SIZE, SIZE, M, ONE, U, IU, JU, DESCU, U, 
00294      $             IU, JU, DESCU, MONE, WORK( PTRR ), 1, 1, DESCR )
00295 *
00296       NORMU = PDLANGE( '1', SIZE, SIZE, WORK( PTRR ), 1, 1, DESCR,
00297      $        WORK( PTRWORK ) )
00298 *
00299       NORMU = NORMU / ULP / SIZE / THRESH
00300       IF( NORMU.GT.1. )
00301      $   RESULT( 2 ) = 1
00302 *
00303 *     Test3. Form identity matrix R  and check norm(VT*VT' - I )
00304 *
00305       CALL PDLASET( 'Full', SIZE, SIZE, ZERO, ONE, WORK( PTRR ), 1, 1,
00306      $              DESCR )
00307       CALL PDGEMM( 'N', 'T', SIZE, SIZE, N, ONE, VT, IVT, JVT, DESCVT,
00308      $             VT, IVT, JVT, DESCVT, MONE, WORK( PTRR ), 
00309      $             1, 1, DESCR )
00310       NORMVT = PDLANGE( '1', SIZE, SIZE, WORK( PTRR ), 1, 1, DESCR,
00311      $         WORK( PTRWORK ) )
00312 *
00313       NORMVT = NORMVT / ULP / SIZE / THRESH
00314       IF( NORMVT.GT.1. )
00315      $   RESULT( 3 ) = 1
00316 *
00317       MTM = MAX( NORMVT, NORMU )*THRESH
00318 *
00319 *     Test 1.
00320 *     Initialize R = diag( S )
00321 *
00322       CALL PDLASET( 'Full', SIZE, SIZE, ZERO, ZERO, WORK( PTRR ), 1, 1,
00323      $              DESCR )
00324 *
00325       DO 10 I = 1, SIZE
00326          CALL PDELSET( WORK( PTRR ), I, I, DESCR, S( I ) )
00327    10 CONTINUE
00328 *
00329 *     Calculate U = U*R
00330 *
00331       DO 20 I = 1, SIZE
00332          PCOL = INDXG2P( I, DESCU( NB_ ), 0, 0, NPCOL )
00333          LOCALCOL = INDXG2L( I, DESCU( NB_ ), 0, 0, NPCOL )
00334          IF( MYCOL.EQ.PCOL ) THEN
00335             CALL DSCAL( MP, S( I ), U( ( LOCALCOL-1 )*DESCU( LLD_ )+1 ),
00336      $                  1 )
00337          END IF
00338    20 CONTINUE
00339 *
00340 *     Calculate A = U*VT - A
00341 *
00342       CALL PDGEMM( 'N', 'N', M, N, SIZE, ONE, U, IU, JU, DESCU, VT,
00343      $             IVT, JVT, DESCVT, MONE, A, IA, JA, DESCA )
00344 *
00345       NORMA = PDLANGE( '1', M, N, A, IA, JA, DESCA, WORK( PTRWORK ) )
00346       THRESHA = NORMAI*MAX( M, N )*ULP*THRESH
00347 *
00348       IF( NORMA.GT.THRESHA )
00349      $   RESULT( 1 ) = 1
00350 *
00351       IF( THRESHA.EQ.0 ) THEN
00352          CHK = 0.0D0
00353       ELSE
00354          CHK = NORMA / THRESHA*THRESH
00355       END IF
00356 *
00357 *     Test 4.
00358 *
00359       DO 30 I = 1, SIZE - 1
00360          FIRST = S( I )
00361          SECOND = S( I+1 )
00362          IF( FIRST.LT.SECOND )
00363      $      RESULT( 4 ) = 1
00364    30 CONTINUE
00365    40 CONTINUE
00366       RETURN
00367       END