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