/***********************************************************************/
/*    LICENSED MATERIALS - PROPERTY OF IBM                             */
/*    "RESTRICTED MATERIALS OF IBM"                                    */
/*                                                                     */
/*    5765-422                                                         */
/*    5765-F84                                                         */
/*    5765-G18                                                         */
/*    (C) COPYRIGHT IBM CORP. 1995, 2003. ALL RIGHTS RESERVED.         */
/*                                                                     */
/*    U.S. GOVERNMENT USERS RESTRICTED RIGHTS - USE, DUPLICATION       */
/*    OR DISCLOSURE RESTRICTED BY GSA ADP SCHEDULE CONTRACT WITH       */
/*    IBM CORP.                                                        */
/***********************************************************************/
#ifdef __linux
  #ifndef min
    #define min(a,b) (((a)<(b))?(a):(b))
  #endif
  #ifndef max
    #define max(a,b) (((a)>(b))?(a):(b))
  #endif
#else
  #include <macros.h>
#endif
#include "parameter.h"
#include <unistd.h>
#include <math.h>
#include <stdio.h>
#include "essl.h"
#include <pessl.h> 


/* external variables defined in the diffusion.c file */

extern int dif_nx, dif_ny, dif_npts, dif_ntemps;
extern double dif_ly_ratio, dif_delta_t, *dif_x, *dif_y;
/***********************************************************************/
/*                                                                     */
/*   dif_nx     contains the number of sine coefficients in x          */
/*   dif_ny     contains the number of sine coefficients in y          */
/*   dif_npts   contains total number of (x,y) solution points to print*/
/*   dif_ntemps contains total number of solution temperatures to print*/
/*   dif_ly_ratio contains the ratio of x dimension to y dimension     */
/*   dif_delta_t contains the temperature step size                    */
/*   dif_x      contains the x coordinates of the printed solution     */
/*   dif_y      contains the y coordinates of the printed solution     */
/*                                                                     */
/***********************************************************************/

void init_diffusion(void);
double init_temp( double x, double y);
double diff_coef( double x, double y);

/* external variables defined in the scale.c file */


extern int sc_icontext, sc_iam, sc_nnodes;
/***********************************************************************/
/*                                                                     */
/*   sc_icontext contains the BLACS context                            */
/*   sc_iam      contains this process's node number                   */
/*   sc_nnodes   contains the total number of nodes                    */
/*                                                                     */
/***********************************************************************/

typedef struct _G_index 
   {
   int num_row_blks;
   int num_col_blks;
   int *g_start_row_blk;
   int *g_start_col_blk;
   int *g_end_row_blk;
   int *g_end_col_blk;
   } G_index;

/* the G_index structure is defined to simplify cases where
   one needs to simultaneously use local and global indices for a
   calculation.  Here the g_end_row_blk is an array, with the size
   num_row_blks, which contains the last global row index in that 
   particular block. The other arrays are correspondingly defined.
*/
void initialize_scale (int n, int *index);
void allocate_rarray( double ***array, int m, int n);
void allocate_carray( dcmplx ***array, int m, int n);
void initialize_rarray( double ***array, int *desc_array, 
                       G_index **index_array, int m, int n, int blk_index);
void initialize_carray( dcmplx ***array, int *desc_array, 
                       G_index **index_array, int m, int n, int blk_index);
        
void clocal_to_rglobal(dcmplx **a, int *a_d, double *a_glob, int lda);
void rlocal_to_rglobal(double **a, int *a_d, double *a_glob, int lda);


void get_diffusion_matrix( double **f, G_index *f_i, int *f_d);
void expand_temp_profile( double **a, G_index *a_i, int *a_d, int lda);

void factor_nodes(void );
int minpower2( int n, int *log2n );


/* define some variable and macros so we can easily use fortran blacs  */
/* routines in c */
int __iblcs1, __iblcs2, __iblcs3, __iblcs4, __iblcs5, __iblcs6, __iblcs7;
double __fblcs1, __fblcs2, __fblcs3, __fblcs4, __fblcs5;

#define blacs_pinfo(x1,x2) (blacs_pinfo(x1, x2))

#define blacs_abort(i1,i2) (__iblcs1=i1, __iblcs2=i2, blacs_abort(&__iblcs1, \
                           &__iblcs2))

#define blacs_get(i1,i2,x1) (__iblcs1=i1, __iblcs2=i2, blacs_get(&__iblcs1, \
                           &__iblcs2,x1))

#define blacs_gridinit(x1,x2,i1,i2) (__iblcs1=i1, __iblcs2=i2, \
                          blacs_gridinit(x1, x2, &__iblcs1, &__iblcs2))

#define blacs_gridinfo(i1,x1,x2,x3,x4) (__iblcs1=i1,             \
                 blacs_gridinfo(&__iblcs1, x1, x2,x3, x4))

#define blacs_barrier(i1,x1) (__iblcs1=i1, blacs_barrier(&__iblcs1, x1))

#define blacs_exit(i1) (__iblcs1=i1, blacs_exit(&__iblcs1))

#define igamx2d(i1,x1,x2,i2,i3,x3,i4,x4,x5,i5,i6,i7)  \
                 ( __iblcs1=i1, __iblcs2=i2,__iblcs3=i3,__iblcs4=i4,  \
                  __iblcs5=i5, __iblcs6=i6,__iblcs7=i7,  \
                  igamx2d(&__iblcs1, x1, x2, &__iblcs2, &__iblcs3, x3, \
                  &__iblcs4, x4, x5, &__iblcs5, &__iblcs6, &__iblcs7))

#define dgebs2d(i1,x1,x2,i2,i3,x3,i4) (__iblcs1=i1, __iblcs2=i2,__iblcs3=i3, \
                  __iblcs4=i4, dgebs2d(&__iblcs1, x1, x2, &__iblcs2, \
                  &__iblcs3, x3,&__iblcs4))

#define dgebr2d(i1,x1,x2,i2,i3,x3,i4,i5,i6) (__iblcs1=i1, __iblcs2=i2,  \
                  __iblcs3=i3, __iblcs4=i4, __iblcs5=i5, __iblcs6=i6, \
                  dgebr2d(&__iblcs1, x1, x2, &__iblcs2, \
                  &__iblcs3, x3,&__iblcs4, &__iblcs5, &__iblcs6))