ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pzgemr2.c
Go to the documentation of this file.
00001 #include "redist.h"
00002 /* $Id: pzgemr2.c,v 1.1.1.1 2000/02/15 18:04:10 susan Exp $
00003  * 
00004  * some functions used by the pzgemr2d routine see file pzgemr.c for more
00005  * documentation.
00006  * 
00007  * Created March 1993 by B. Tourancheau (See sccs for modifications). */
00008 #define static2 static
00009 #if defined(Add_) || defined(f77IsF2C)
00010 #define fortran_mr2d pzgemr2do_
00011 #define fortran_mr2dnew pzgemr2d_
00012 #elif defined(UpCase)
00013 #define fortran_mr2dnew PZGEMR2D
00014 #define fortran_mr2d PZGEMR2DO
00015 #define zcopy_ ZCOPY
00016 #define zlacpy_ ZLACPY
00017 #else
00018 #define fortran_mr2d pzgemr2do
00019 #define fortran_mr2dnew pzgemr2d
00020 #define zcopy_ zcopy
00021 #define zlacpy_ zlacpy
00022 #endif
00023 #define Clacpy Czgelacpy
00024 void  Clacpy();
00025 typedef struct {
00026   double r, i;
00027 }     dcomplex;
00028 typedef struct {
00029   int   desctype;
00030   int   ctxt;
00031   int   m;
00032   int   n;
00033   int   nbrow;
00034   int   nbcol;
00035   int   sprow;
00036   int   spcol;
00037   int   lda;
00038 }     MDESC;
00039 #define BLOCK_CYCLIC_2D 1
00040 typedef struct {
00041   int   lstart;
00042   int   len;
00043 }     IDESC;
00044 #define SHIFT(row,sprow,nbrow) ((row)-(sprow)+ ((row) >= (sprow) ? 0 : (nbrow)))
00045 #define max(A,B) ((A)>(B)?(A):(B))
00046 #define min(A,B) ((A)>(B)?(B):(A))
00047 #define DIVUP(a,b) ( ((a)-1) /(b)+1)
00048 #define ROUNDUP(a,b) (DIVUP(a,b)*(b))
00049 #ifdef MALLOCDEBUG
00050 #define malloc mymalloc
00051 #define free myfree
00052 #define realloc myrealloc
00053 #endif
00054 /* Cblacs */
00055 extern void Cblacs_pcoord();
00056 extern int Cblacs_pnum();
00057 extern void Csetpvmtids();
00058 extern void Cblacs_get();
00059 extern void Cblacs_pinfo();
00060 extern void Cblacs_gridinfo();
00061 extern void Cblacs_gridinit();
00062 extern void Cblacs_exit();
00063 extern void Cblacs_gridexit();
00064 extern void Cblacs_setup();
00065 extern void Cigebs2d();
00066 extern void Cigebr2d();
00067 extern void Cigesd2d();
00068 extern void Cigerv2d();
00069 extern void Cigsum2d();
00070 extern void Cigamn2d();
00071 extern void Cigamx2d();
00072 extern void Czgesd2d();
00073 extern void Czgerv2d();
00074 /* lapack */
00075 void  zlacpy_();
00076 /* aux fonctions */
00077 extern int localindice();
00078 extern void *mr2d_malloc();
00079 extern int ppcm();
00080 extern int localsize();
00081 extern int memoryblocksize();
00082 extern int changeorigin();
00083 extern void paramcheck();
00084 /* tools and others function */
00085 #define scanD0 zgescanD0
00086 #define dispmat zgedispmat
00087 #define setmemory zgesetmemory
00088 #define freememory zgefreememory
00089 #define scan_intervals zgescan_intervals
00090 extern void scanD0();
00091 extern void dispmat();
00092 extern void setmemory();
00093 extern void freememory();
00094 extern int scan_intervals();
00095 extern void Cpzgemr2do();
00096 extern void Cpzgemr2d();
00097 /* some defines for Cpzgemr2do */
00098 #define SENDBUFF 0
00099 #define RECVBUFF 1
00100 #define SIZEBUFF 2
00101 #if 0
00102 #define DEBUG
00103 #endif
00104 #ifndef DEBUG
00105 #define NDEBUG
00106 #endif
00107 #include <stdio.h>
00108 #include <stdlib.h>
00109 #include <string.h>
00110 #include <assert.h>
00111 #include <ctype.h>
00112 /* Created March 1993 by B. Tourancheau (See sccs for modifications). */
00113 /************************************************************************/
00114 /* Set the memory space with the malloc function */
00115 void
00116 setmemory(adpointer, blocksize)
00117   dcomplex **adpointer;
00118   int   blocksize;
00119 {
00120   assert(blocksize >= 0);
00121   if (blocksize == 0) {
00122     *adpointer = NULL;
00123     return;
00124   }
00125   *adpointer = (dcomplex *) mr2d_malloc(
00126                                         blocksize * sizeof(dcomplex));
00127 }
00128 /******************************************************************/
00129 /* Free the memory space after the malloc */
00130 void
00131 freememory(ptrtobefreed)
00132   dcomplex *ptrtobefreed;
00133 {
00134   if (ptrtobefreed == NULL)
00135     return;
00136   free((char *) ptrtobefreed);
00137 }
00138 /* extern functions for intersect() extern zcopy_(); */
00139 /* scan_intervals: scans two distributions in one dimension, and compute the
00140  * intersections on the local processor. result must be long enough to
00141  * contains the result that are stocked in IDESC structure, the function
00142  * returns the number of intersections found */
00143 int 
00144 scan_intervals(type, ja, jb, n, ma, mb, q0, q1, col0, col1,
00145                result)
00146   char  type;
00147   int   ja, jb, n, q0, q1, col0, col1;
00148   MDESC *ma, *mb;
00149   IDESC *result;
00150 {
00151   int   offset, j0, j1, templatewidth0, templatewidth1, nbcol0, nbcol1;
00152   int   l;      /* local indice on the beginning of the interval */
00153   assert(type == 'c' || type == 'r');
00154   nbcol0 = (type == 'c' ? ma->nbcol : ma->nbrow);
00155   nbcol1 = (type == 'c' ? mb->nbcol : mb->nbrow);
00156   templatewidth0 = q0 * nbcol0;
00157   templatewidth1 = q1 * nbcol1;
00158   {
00159     int   sp0 = (type == 'c' ? ma->spcol : ma->sprow);
00160     int   sp1 = (type == 'c' ? mb->spcol : mb->sprow);
00161     j0 = SHIFT(col0, sp0, q0) * nbcol0 - ja;
00162     j1 = SHIFT(col1, sp1, q1) * nbcol1 - jb;
00163   }
00164   offset = 0;
00165   l = 0;
00166   /* a small check to verify that the submatrix begin inside the first block
00167    * of the original matrix, this done by a sort of coordinate change at the
00168    * beginning of the Cpzgemr2d */
00169   assert(j0 + nbcol0 > 0);
00170   assert(j1 + nbcol1 > 0);
00171   while ((j0 < n) && (j1 < n)) {
00172     int   end0, end1;
00173     int   start, end;
00174     end0 = j0 + nbcol0;
00175     end1 = j1 + nbcol1;
00176     if (end0 <= j1) {
00177       j0 += templatewidth0;
00178       l += nbcol0;
00179       continue;
00180     }
00181     if (end1 <= j0) {
00182       j1 += templatewidth1;
00183       continue;
00184     }
00185     /* compute the raw intersection */
00186     start = max(j0, j1);
00187     start = max(start, 0);
00188     /* the start is correct now, update the corresponding fields */
00189     result[offset].lstart = l + start - j0;
00190     end = min(end0, end1);
00191     if (end0 == end) {
00192       j0 += templatewidth0;
00193       l += nbcol0;
00194     }
00195     if (end1 == end)
00196       j1 += templatewidth1;
00197     /* throw the limit if they go out of the matrix */
00198     end = min(end, n);
00199     assert(end > start);
00200     /* it is a bit tricky to see why the length is always positive after all
00201      * this min and max, first we have the property that every interval
00202      * considered is at least partly into the submatrix, second we arrive
00203      * here only if the raw intersection is non-void, if we remove a limit
00204      * that means the corresponding frontier is in both intervals which
00205      * proove the final interval is non-void, clear ?? */
00206     result[offset].len = end - start;
00207     offset += 1;
00208   }     /* while */
00209   return offset;
00210 }