ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
PB_Cplasca2.c
Go to the documentation of this file.
00001 /* ---------------------------------------------------------------------
00002 *
00003 *  -- PBLAS auxiliary routine (version 2.0) --
00004 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00005 *     and University of California, Berkeley.
00006 *     April 1, 1998
00007 *
00008 *  ---------------------------------------------------------------------
00009 */
00010 /*
00011 *  Include files
00012 */
00013 #include "../pblas.h"
00014 #include "../PBpblas.h"
00015 #include "../PBtools.h"
00016 #include "../PBblacs.h"
00017 #include "../PBblas.h"
00018 
00019 #ifdef __STDC__
00020 void PB_Cplasca2( PBTYP_T * TYPE, char * UPLO, char * CONJUG, int M,
00021                   int N, char * ALPHA, char * A, int IA, int JA,
00022                   int * DESCA )
00023 #else
00024 void PB_Cplasca2( TYPE, UPLO, CONJUG, M, N, ALPHA, A, IA, JA, DESCA )
00025 /*
00026 *  .. Scalar Arguments ..
00027 */
00028    char           * CONJUG, * UPLO;
00029    int            IA, JA, M, N;
00030    char           * ALPHA;
00031    PBTYP_T        * TYPE;
00032 /*
00033 *  .. Array Arguments ..
00034 */
00035    int            * DESCA;
00036    char           * A;
00037 #endif
00038 {
00039 /*
00040 *  .. Local Scalars ..
00041 */
00042    char           UploA, herm;
00043    int            Acol, Arow, Aii, iimax, ilow, imbloc, Aimb1, inbloc, Ainb1,
00044                   Aoffi, GoEast, GoSouth, ioffd, iupp, izero=0, Ajj, jjmax,
00045                   Aoffj, joffd, lcmt, lcmt00, Ald, lmbloc, lnbloc, low, lower,
00046                   m1, Amb, mbloc, mblkd, mblks, Amp, Arcol, Arrow, mycol, myrow,
00047                   n1, Anb, nbloc, nblkd, nblks, npcol, nprow, Anq, pmb, qnb,
00048                   size, tmp1, upp, upper;
00049    TZSCAL_T       scal;
00050 /* ..
00051 *  .. Executable Statements ..
00052 *
00053 */
00054 /*
00055 *  Quick return if possible
00056 */
00057    if( ( M <= 0 ) || ( N <= 0 ) ) return;
00058 /*
00059 *  Retrieve process grid information
00060 */
00061    Cblacs_gridinfo( DESCA[CTXT_], &nprow, &npcol, &myrow, &mycol );
00062 /*
00063 *  Retrieve sub( A )'s local information: Aii, Ajj, Arow, Acol ...
00064 */
00065    PB_Cainfog2l( M, N, IA, JA, DESCA, nprow, npcol, myrow, mycol, &Aimb1,
00066                  &Ainb1, &Amp, &Anq, &Aii, &Ajj, &Arow, &Acol, &Arrow, &Arcol );
00067 /*
00068 *  Quick return if I don't own any of sub( A ).
00069 */
00070    if( ( Amp <= 0 ) || ( Anq <= 0 ) ) return;
00071 /*
00072 *  Initialize lcmt00, mblks, nblks, imbloc, inbloc, lmbloc, lnbloc, ilow, low,
00073 *  iupp, and upp.
00074 */
00075    Amb   = DESCA[MB_ ]; Anb   = DESCA[NB_ ]; Ald   = DESCA[LLD_];
00076    PB_Cbinfo( 0, Amp, Anq, Aimb1, Ainb1, Amb, Anb, Arrow, Arcol, &lcmt00,
00077               &mblks, &nblks, &imbloc, &inbloc, &lmbloc, &lnbloc, &ilow, &low,
00078               &iupp, &upp );
00079    iimax = ( Aoffi = Aii - 1 ) + ( m1 = Amp );
00080    jjmax = ( Aoffj = Ajj - 1 ) + ( n1 = Anq );
00081    pmb = ( ( ( Arow < 0 ) || ( nprow == 1 ) ) ? Amb : nprow * Amb );
00082    qnb = ( ( ( Acol < 0 ) || ( npcol == 1 ) ) ? Anb : npcol * Anb );
00083 
00084    UploA = Mupcase( UPLO[0] );
00085    upper = ( UploA != CLOWER );
00086    lower = ( UploA != CUPPER );
00087    herm  = ( UploA == CALL ? CNOCONJG : Mupcase( CONJUG[0] ) );
00088 
00089    size  = TYPE->size;
00090    scal  = ( herm == CCONJG ? TYPE->Fhescal : TYPE->Ftzscal );
00091 /*
00092 *  Handle separately the first row and/or column of the LCM table. Update the
00093 *  LCM value of the curent block lcmt00, as well as the number of rows and
00094 *  columns mblks and nblks remaining in the LCM table.
00095 */
00096    GoSouth = ( lcmt00 > iupp );
00097    GoEast  = ( lcmt00 < ilow );
00098 /*
00099 *  Go through the table looking for blocks owning diagonal entries.
00100 */
00101    if( ( !( GoSouth ) ) && ( !( GoEast ) ) )
00102    {
00103 /*
00104 *  The upper left block owns diagonal entries lcmt00 >= ilow && lcmt00 <= iupp
00105 */
00106       scal( C2F_CHAR( UPLO ), &imbloc, &inbloc, &lcmt00, ALPHA,
00107             Mptr( A, Aii, Ajj, Ald, size ), &Ald );
00108 /*
00109 *  Decide whether one should go south or east in the table: Go east if
00110 *  the block below the current one only owns lower entries. If this block,
00111 *  however, owns diagonals, then go south.
00112 */
00113       GoSouth = !( GoEast = ( ( lcmt00 - ( iupp - upp + pmb ) ) < ilow ) );
00114 
00115       if( GoSouth )
00116       {
00117 /*
00118 *  When the upper triangular part of sub( A ) should be scaled and one is
00119 *  planning to go south in the table, it is neccessary to take care of the
00120 *  remaining columns of these imbloc rows immediately.
00121 */
00122          if( upper && ( Anq > inbloc ) )
00123          {
00124             tmp1 = Anq - inbloc;
00125             scal( C2F_CHAR( ALL ), &imbloc, &tmp1, &izero, ALPHA,
00126                   Mptr( A, Aii, Ajj+inbloc, Ald, size ), &Ald );
00127          }
00128          Aii += imbloc;
00129          m1  -= imbloc;
00130       }
00131       else
00132       {
00133 /*
00134 *  When the lower triangular part of sub( A ) should be scaled and one is
00135 *  planning to go east in the table, it is neccessary to take care of the
00136 *  remaining rows of these inbloc columns immediately.
00137 */
00138          if( lower && ( Amp > imbloc ) )
00139          {
00140             tmp1 = Amp - imbloc;
00141             scal( C2F_CHAR( ALL ), &tmp1, &inbloc, &izero, ALPHA,
00142                   Mptr( A, Aii+imbloc, Ajj, Ald, size ), &Ald );
00143          }
00144          Ajj += inbloc;
00145          n1  -= inbloc;
00146       }
00147    }
00148 
00149    if( GoSouth )
00150    {
00151 /*
00152 *  Go one step south in the LCM table. Adjust the current LCM value as well as
00153 *  the local row index in A.
00154 */
00155       lcmt00 -= ( iupp - upp + pmb ); mblks--; Aoffi += imbloc;
00156 /*
00157 *  While there are blocks remaining that own upper entries, keep going south.
00158 *  Adjust the current LCM value as well as the local row index in A.
00159 */
00160       while( ( mblks > 0 ) && ( lcmt00 > upp ) )
00161       { lcmt00 -= pmb; mblks--; Aoffi += Amb; }
00162 /*
00163 *  Scale the upper triangular part of sub( A ) we just skipped when necessary.
00164 */
00165       tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
00166       if( upper && ( tmp1 > 0 ) )
00167       {
00168          scal( C2F_CHAR( ALL ), &tmp1, &n1, &izero, ALPHA,
00169                Mptr( A, Aii, Aoffj+1, Ald, size ), &Ald );
00170          Aii += tmp1;
00171          m1  -= tmp1;
00172       }
00173 /*
00174 *  Return if no more row in the LCM table.
00175 */
00176       if( mblks <= 0 ) return;
00177 /*
00178 *  lcmt00 <= upp. The current block owns either diagonals or lower entries.
00179 *  Save the current position in the LCM table. After this column has been
00180 *  completely taken care of, re-start from this row and the next column of
00181 *  the LCM table.
00182 */
00183       lcmt  = lcmt00; mblkd = mblks; ioffd = Aoffi;
00184 
00185       mbloc = Amb;
00186       while( ( mblkd > 0 ) && ( lcmt >= ilow ) )
00187       {
00188 /*
00189 *  A block owning diagonals lcmt00 >= ilow && lcmt00 <= upp has been found.
00190 */
00191          if( mblkd == 1 ) mbloc = lmbloc;
00192          scal( C2F_CHAR( UPLO ), &mbloc, &inbloc, &lcmt, ALPHA,
00193                Mptr( A, ioffd+1, Aoffj+1, Ald, size ), &Ald );
00194          lcmt00  = lcmt;
00195          lcmt   -= pmb;
00196          mblks   = mblkd;
00197          mblkd--;
00198          Aoffi   = ioffd;
00199          ioffd  += mbloc;
00200       }
00201 /*
00202 *  Scale the lower triangular part of sub( A ) when necessary.
00203 */
00204       tmp1 = m1 - ioffd + Aii - 1;
00205       if( lower && ( tmp1 > 0 ) )
00206          scal( C2F_CHAR( ALL ), &tmp1, &inbloc, &izero, ALPHA,
00207                Mptr( A, ioffd+1, Aoffj+1, Ald, size ), &Ald );
00208 
00209       tmp1    = Aoffi - Aii + 1;
00210       m1     -= tmp1;
00211       n1     -= inbloc;
00212       lcmt00 += low - ilow + qnb;
00213       nblks--;
00214       Aoffj  += inbloc;
00215 /*
00216 *  When the upper triangular part of sub( A ) should be scaled, take care of the
00217 *  n1 remaining columns of these tmp1 rows immediately.
00218 */
00219       if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) )
00220          scal( C2F_CHAR( ALL ), &tmp1, &n1, &izero, ALPHA,
00221                Mptr( A, Aii, Aoffj+1, Ald, size ), &Ald );
00222       Aii = Aoffi + 1;
00223       Ajj = Aoffj + 1;
00224    }
00225    else if( GoEast )
00226    {
00227 /*
00228 *  Go one step east in the LCM table. Adjust the current LCM value as well as
00229 *  the local column index in A.
00230 */
00231       lcmt00 += low - ilow + qnb; nblks--; Aoffj += inbloc;
00232 /*
00233 *  While there are blocks remaining that own lower entries, keep going east.
00234 *  Adjust the current LCM value as well as the local column index in A.
00235 */
00236       while( ( nblks > 0 ) && ( lcmt00 < low ) )
00237       { lcmt00 += qnb; nblks--; Aoffj += Anb; }
00238 /*
00239 *  Scale the lower triangular part of sub( A ) we just skipped when necessary.
00240 */
00241       tmp1 = MIN( Aoffj, jjmax ) - Ajj + 1;
00242       if( lower && ( tmp1 > 0 ) )
00243       {
00244          scal( C2F_CHAR( ALL ), &m1, &tmp1, &izero, ALPHA,
00245                Mptr( A, Aii, Ajj, Ald, size ), &Ald );
00246          Ajj += tmp1;
00247          n1  -= tmp1;
00248       }
00249 /*
00250 *  Return if no more column in the LCM table.
00251 */
00252       if( nblks <= 0 ) return;
00253 /*
00254 *  lcmt00 >= low. The current block owns either diagonals or upper entries.
00255 *  Save the current position in the LCM table. After this row has been
00256 *  completely taken care of, re-start from this column and the next row of
00257 *  the LCM table.
00258 */
00259       lcmt  = lcmt00; nblkd = nblks; joffd = Aoffj;
00260 
00261       nbloc = Anb;
00262       while( ( nblkd > 0 ) && ( lcmt <= iupp ) )
00263       {
00264 /*
00265 *  A block owning diagonals lcmt00 >= low && lcmt00 <= iupp has been found.
00266 */
00267          if( nblkd == 1 ) nbloc = lnbloc;
00268          scal( C2F_CHAR( UPLO ), &imbloc, &nbloc, &lcmt, ALPHA,
00269                Mptr( A, Aii, joffd+1, Ald, size ), &Ald );
00270          lcmt00  = lcmt;
00271          lcmt   += qnb;
00272          nblks   = nblkd;
00273          nblkd--;
00274          Aoffj   = joffd;
00275          joffd  += nbloc;
00276       }
00277 /*
00278 *  Scale the upper triangular part of sub( A ) when necessary.
00279 */
00280       tmp1 = n1 - joffd + Ajj - 1;
00281       if( upper && ( tmp1 > 0 ) )
00282          scal( C2F_CHAR( ALL ), &imbloc, &tmp1, &izero, ALPHA,
00283                Mptr( A, Aii, joffd+1, Ald, size ), &Ald );
00284 
00285       tmp1    = Aoffj - Ajj + 1;
00286       m1     -= imbloc;
00287       n1     -= tmp1;
00288       lcmt00 -= ( iupp - upp + pmb );
00289       mblks--;
00290       Aoffi  += imbloc;
00291 /*
00292 *  When the lower triangular part of sub( A ) should be scaled, take care of the
00293 *  m1 remaining rows of these tmp1 columns immediately.
00294 */
00295       if( lower && ( m1 > 0 ) && ( tmp1 > 0 ) )
00296          scal( C2F_CHAR( ALL ), &m1, &tmp1, &izero, ALPHA,
00297                Mptr( A, Aoffi+1, Ajj, Ald, size ), &Ald );
00298       Aii = Aoffi + 1;
00299       Ajj = Aoffj + 1;
00300    }
00301 /*
00302 *  Loop over the remaining columns of the LCM table.
00303 */
00304    nbloc = Anb;
00305    while( nblks > 0 )
00306    {
00307       if( nblks == 1 ) nbloc = lnbloc;
00308 /*
00309 *  While there are blocks remaining that own upper entries, keep going south.
00310 *  Adjust the current LCM value as well as the local row index in A.
00311 */
00312       while( ( mblks > 0 ) && ( lcmt00 > upp ) )
00313       { lcmt00 -= pmb; mblks--; Aoffi  += Amb; }
00314 /*
00315 *  Scale the upper triangular part of sub( A ) we just skipped when necessary.
00316 */
00317       tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
00318       if( upper && ( tmp1 > 0 ) )
00319       {
00320          scal( C2F_CHAR( ALL ), &tmp1, &n1, &izero, ALPHA,
00321                Mptr( A, Aii, Aoffj+1, Ald, size ), &Ald );
00322          Aii += tmp1;
00323          m1  -= tmp1;
00324       }
00325 /*
00326 *  Return if no more row in the LCM table.
00327 */
00328       if( mblks <= 0 ) return;
00329 /*
00330 *  lcmt00 <= upp. The current block owns either diagonals or lower entries.
00331 *  Save the current position in the LCM table. After this column has been
00332 *  completely taken care of, re-start from this row and the next column of
00333 *  the LCM table.
00334 */
00335       lcmt  = lcmt00; mblkd = mblks; ioffd = Aoffi;
00336 
00337       mbloc = Amb;
00338       while( ( mblkd > 0 ) && ( lcmt >= low ) )
00339       {
00340 /*
00341 *  A block owning diagonals lcmt00 >= low && lcmt00 <= upp has been found.
00342 */
00343          if( mblkd == 1 ) mbloc = lmbloc;
00344          scal( C2F_CHAR( UPLO ), &mbloc, &nbloc, &lcmt, ALPHA,
00345                Mptr( A, ioffd+1, Aoffj+1, Ald, size ), &Ald );
00346          lcmt00  = lcmt;
00347          lcmt   -= pmb;
00348          mblks   = mblkd;
00349          mblkd--;
00350          Aoffi   = ioffd;
00351          ioffd  += mbloc;
00352       }
00353 /*
00354 *  Scale the lower triangular part of sub( A ) when necessary.
00355 */
00356       tmp1 = m1 - ioffd + Aii - 1;
00357       if( lower && ( tmp1 > 0 ) )
00358          scal( C2F_CHAR( ALL ), &tmp1, &nbloc, &izero, ALPHA,
00359                Mptr( A, ioffd+1, Aoffj+1, Ald, size ), &Ald );
00360 
00361       tmp1    = MIN( Aoffi, iimax ) - Aii + 1;
00362       m1     -= tmp1;
00363       n1     -= nbloc;
00364       lcmt00 += qnb;
00365       nblks--;
00366       Aoffj  += nbloc;
00367 /*
00368 *  When the upper triangular part of sub( A ) should be scaled, take care of the
00369 *  n1 remaining columns of these tmp1 rows immediately.
00370 */
00371       if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) )
00372          scal( C2F_CHAR( ALL ), &tmp1, &n1, &izero, ALPHA,
00373                Mptr( A, Aii, Aoffj+1, Ald, size ), &Ald );
00374       Aii = Aoffi + 1;
00375       Ajj = Aoffj + 1;
00376    }
00377 /*
00378 *  End of PB_Cplasca2
00379 */
00380 }