← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

Viewing changes to src/tools/ga-5-1/global/testing/elempatch.c

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools/ga-5-1/global/testing/elempatch.c
1
 
#if HAVE_CONFIG_H
2
 
#   include "config.h"
3
 
#endif
4
 
 
5
 
/**\
6
 
File: elempatch.c
7
 
Purpose: test the interfaces:
8
 
 
9
 
    GA_Abs_value_patch(g_a)
10
 
    GA_Add_constant_patch(g_a, alpha)
11
 
    GA_Recip_patch_patch(g_a)
12
 
    GA_Elem_multiply_patch(g_a, alo, ahi, g_b, blo, bhi, g_c, clo, chi)
13
 
    GA_Elem_divide_patch(g_a, alo, ahi, g_b, blo, bhi, g_c, clo, chi)
14
 
    GA_Elem_maximum_patch(g_a, alo, ahi, g_b, blo, bhi, g_c, clo, ch
15
 
    GA_Elem_minimum_patch(g_a, alo, ahi, g_b, blo, bhi, g_c, clo, chi)
16
 
        
17
 
        that are for TAO/Global Array Project
18
 
 
19
 
Author:
20
 
 
21
 
Limin Zhang, Ph.D.
22
 
Mathematics Department
23
 
Columbia Basin College
24
 
Pasco, WA 99301
25
 
 
26
 
Mentor:
27
 
 
28
 
Jarek Nieplocha
29
 
Pacific Northwest National Laboratory
30
 
 
31
 
Date: Jauary 30, 2002
32
 
Revised on February 26, 2002.
33
 
 
34
 
\**/
35
 
 
36
 
#if HAVE_MATH_H
37
 
#   include <math.h>
38
 
#endif
39
 
#if HAVE_STDLIB_H
40
 
#   include <stdlib.h>
41
 
#endif
42
 
#if HAVE_STDIO_H
43
 
#   include <stdio.h>
44
 
#endif
45
 
 
46
 
#include "ga.h"
47
 
#include "macdecls.h"
48
 
#include "mp3.h"
49
 
#include "../src/globalp.h"
50
 
#include "ga-papi.h"
51
 
 
52
 
#define BLOCK_CYCLIC 0
53
 
 
54
 
#if BLOCK_CYCLIC
55
 
#define USE_SCALAPACK 1
56
 
#endif
57
 
 
58
 
#ifndef GA_HALF_MAX_INT 
59
 
#define GA_HALF_MAX_INT ((((int)1) << ((int)(8*sizeof(int))-2)) - 1)
60
 
#endif
61
 
 
62
 
#ifndef GA_INFINITY_I
63
 
#define GA_INFINITY_I (GA_HALF_MAX_INT + GA_HALF_MAX_INT + 1)
64
 
/* 
65
 
  Original value below.
66
 
  Seemed too small arbitrarily.
67
 
#define GA_INFINITY_I 100000
68
 
*/
69
 
#endif
70
 
 
71
 
#ifndef GA_NEGATIVE_INFINITY_I
72
 
#define GA_NEGATIVE_INFINITY_I (- GA_INFINITY_I)
73
 
 
74
 
 
75
 
/* 
76
 
  Original value below. 
77
 
  Seemed too small arbitrarily.
78
 
#define GA_NEGATIVE_INFINITY_I -100000
79
 
*/
80
 
#endif
81
 
 
82
 
#ifndef GA_HALF_MAX_LONG
83
 
#define GA_HALF_MAX_LONG ((((long)1) << ((int)(8*sizeof(long))-2)) - 1)
84
 
#endif
85
 
 
86
 
#ifndef GA_INFINITY_L
87
 
#define GA_INFINITY_L (GA_HALF_MAX_LONG + GA_HALF_MAX_LONG + 1)
88
 
/* Original value was
89
 
#define GA_INFINITY_L 100000
90
 
*/
91
 
#endif
92
 
 
93
 
#ifndef GA_NEGATIVE_INFINITY_L
94
 
#define GA_NEGATIVE_INFINITY_L (- GA_INFINITY_L)
95
 
#endif
96
 
/* 
97
 
  Original value was:
98
 
#define GA_NEGATIVE_INFINITY_L -100000
99
 
*/
100
 
 
101
 
/* 
102
 
  Modified by Doug Baxter 01/24/04 to distinguish between
103
 
  Double inifinity and float infinity.
104
 
#ifndef GA_INFINITY
105
 
#define GA_INFINITY 1.0e20
106
 
#endif
107
 
 
108
 
#ifndef GA_NEGATIVE_INFINITY
109
 
#define GA_NEGATIVE_INFINITY -1.0e20
110
 
#endif
111
 
*/
112
 
#ifndef GA_INFINITY_F
113
 
#define GA_INFINITY_F 1.0e37
114
 
#endif
115
 
/*
116
 
  Original value below.
117
 
#define GA_INFINITY_F 1.0e20
118
 
*/
119
 
#ifndef GA_NEGATIVE_INFINITY_F
120
 
#define GA_NEGATIVE_INFINITY_F -1.0e37
121
 
#endif
122
 
/*
123
 
  Original value below.
124
 
#define GA_NEGATIVE_INFINITY_F -1.0e20
125
 
*/
126
 
#ifndef GA_INFINITY_D
127
 
#define GA_INFINITY_D 1.0e307
128
 
#endif
129
 
/*
130
 
  Original value below.
131
 
#define GA_INFINITY_D 1.0e20
132
 
*/
133
 
#ifndef GA_NEGATIVE_INFINITY_D
134
 
#define GA_NEGATIVE_INFINITY_D -1.0e307
135
 
#endif
136
 
 
137
 
 
138
 
# define THRESH 1e-5
139
 
#define MISMATCHED(x,y) GA_ABS((x)-(y))>=THRESH
140
 
 
141
 
#define N 100
142
 
#define BLOCK_SIZE 20
143
 
#define OP_ELEM_MULT 0
144
 
#define OP_ELEM_DIV 1
145
 
#define OP_ELEM_MAX 2
146
 
#define OP_ELEM_MIN 3
147
 
#define OP_ABS 4
148
 
#define OP_ADD_CONST 5
149
 
#define OP_RECIP 6
150
 
#define OP_STEP_MAX 7
151
 
#define OP_STEP_BOUND_INFO 8
152
 
#define MY_TYPE 2002
153
 
 
154
 
Integer _ga_lo[MAXDIM], _ga_hi[MAXDIM], _ga_work[MAXDIM];
155
 
#  define COPYINDEX_C2F(carr, farr, n){\
156
 
   int i; for(i=0; i< (n); i++)(farr)[n-i-1]=(Integer)(carr)[i]+1;}
157
 
 
158
 
void nga_vfill_patch(Integer *g_a, Integer *lo, Integer *hi);
159
 
void nga_pnfill_patch(Integer *g_a, Integer *lo, Integer *hi);
160
 
 
161
 
void NGA_Vfill_patch(int g_a, int lo[], int hi[])
162
 
{
163
 
    Integer a=(Integer)g_a;
164
 
    Integer ndim = pnga_ndim(a);
165
 
    COPYINDEX_C2F(lo,_ga_lo, ndim);
166
 
    COPYINDEX_C2F(hi,_ga_hi, ndim);
167
 
 
168
 
    nga_vfill_patch(&a, _ga_lo, _ga_hi);
169
 
}
170
 
 
171
 
 
172
 
void NGA_Pnfill_patch(int g_a, int lo[], int hi[])
173
 
{
174
 
    Integer a=(Integer)g_a;
175
 
    Integer ndim = pnga_ndim(a);
176
 
    COPYINDEX_C2F(lo,_ga_lo, ndim);
177
 
    COPYINDEX_C2F(hi,_ga_hi, ndim);
178
 
 
179
 
    nga_pnfill_patch(&a, _ga_lo, _ga_hi);
180
 
}
181
 
 
182
 
int
183
 
ifun (int k)
184
 
{
185
 
  int result;
186
 
  result = -k - 1;
187
 
  result = -2;
188
 
  return result;
189
 
}
190
 
 
191
 
int
192
 
ifun2 (int k)
193
 
{
194
 
  int result;
195
 
  result = k + 1;
196
 
  result = -3;
197
 
  return result;
198
 
}
199
 
 
200
 
void
201
 
fill_func (int nelem, int type, void *buf)
202
 
{
203
 
  int i;
204
 
 
205
 
 
206
 
  switch (type)
207
 
    {
208
 
    case C_FLOAT:
209
 
      for (i = 0; i < nelem; i++)
210
 
    ((float *) buf)[i] = (float) ifun (i);
211
 
      break;
212
 
    case C_LONG:
213
 
      for (i = 0; i < nelem; i++)
214
 
    ((long *) buf)[i] = (long) ifun (i);
215
 
      break;
216
 
    case C_DBL:
217
 
      for (i = 0; i < nelem; i++)
218
 
    ((double *) buf)[i] = (double) ifun (i);
219
 
      break;
220
 
    case C_DCPL:
221
 
      for (i = 0; i < 2 * nelem; i++)
222
 
    ((double *) buf)[i] = (double) ifun (i);
223
 
      break;
224
 
    case C_SCPL:
225
 
      for (i = 0; i < 2 * nelem; i++)
226
 
    ((float *) buf)[i] = (float) ifun (i);
227
 
      break;
228
 
    case C_INT:
229
 
      for (i = 0; i < nelem; i++)
230
 
    ((int *) buf)[i] = ifun (i);
231
 
      break;
232
 
    default:
233
 
      GA_Error (" wrong data type ", type);
234
 
 
235
 
    }
236
 
}
237
 
 
238
 
void
239
 
fill_func2 (int nelem, int type, void *buf)
240
 
{
241
 
  /* int i,size=MA_sizeof(MT_CHAR,type,1);*/
242
 
 
243
 
  int i;
244
 
 
245
 
  switch (type)
246
 
    {
247
 
    case C_FLOAT:
248
 
      for (i = 0; i < nelem; i++)
249
 
    ((float *) buf)[i] = (float) ifun2 (i);
250
 
      break;
251
 
    case C_LONG:
252
 
      for (i = 0; i < nelem; i++)
253
 
    ((long *) buf)[i] = (long) ifun2 (i);
254
 
      break;
255
 
    case C_DBL:
256
 
      for (i = 0; i < nelem; i++)
257
 
    ((double *) buf)[i] = (double) ifun2 (i);
258
 
      break;
259
 
    case C_DCPL:
260
 
      for (i = 0; i < 2 * nelem; i++)
261
 
    ((double *) buf)[i] = (double) ifun2 (i);
262
 
      break;
263
 
    case C_SCPL:
264
 
      for (i = 0; i < 2 * nelem; i++)
265
 
    ((float *) buf)[i] = (float) ifun2 (i);
266
 
      break;
267
 
    case C_INT:
268
 
      for (i = 0; i < nelem; i++)
269
 
    ((int *) buf)[i] = ifun2 (i);
270
 
      break;
271
 
    default:
272
 
      GA_Error (" wrong data type ", type);
273
 
 
274
 
    }
275
 
}
276
 
 
277
 
void
278
 
fill_func3 (int nelem, int type, void *buf)
279
 
/*taking the absolute of the ifun() */
280
 
{
281
 
/*int i,size=MA_sizeof(MT_CHAR,type,1);*/
282
 
 
283
 
  int i;
284
 
 
285
 
  switch (type)
286
 
    {
287
 
    case C_FLOAT:
288
 
      for (i = 0; i < nelem; i++)
289
 
    ((float *) buf)[i] = (float) GA_ABS (ifun (i));
290
 
      break;
291
 
    case C_LONG:
292
 
      for (i = 0; i < nelem; i++)
293
 
    ((long *) buf)[i] = (long) GA_ABS (ifun (i));
294
 
      break;
295
 
    case C_DBL:
296
 
      for (i = 0; i < nelem; i++)
297
 
    ((double *) buf)[i] = (double) GA_ABS (ifun (i));
298
 
      break;
299
 
    case C_DCPL:
300
 
      for (i = 0; i < 2 * nelem - 1; i = i + 2)
301
 
    {
302
 
      ((double *) buf)[i] =
303
 
        sqrt ((double)
304
 
          (ifun (i) * ifun (i) + ifun (i + 1) * ifun (i + 1)));
305
 
      ((double *) buf)[i + 1] = 0.0;
306
 
    }
307
 
      break;
308
 
    case C_SCPL:
309
 
      for (i = 0; i < 2 * nelem - 1; i = i + 2)
310
 
    {
311
 
      ((float *) buf)[i] =
312
 
        sqrt ((float)
313
 
          (ifun (i) * ifun (i) + ifun (i + 1) * ifun (i + 1)));
314
 
      ((float *) buf)[i + 1] = 0.0;
315
 
    }
316
 
      break;
317
 
    case C_INT:
318
 
      for (i = 0; i < nelem; i++)
319
 
    ((int *) buf)[i] = GA_ABS (ifun (i));
320
 
      break;
321
 
    default:
322
 
      GA_Error (" wrong data type ", type);
323
 
 
324
 
    }
325
 
}
326
 
 
327
 
 
328
 
 
329
 
 
330
 
 
331
 
 
332
 
int
333
 
test_fun (int type, int dim, int OP)
334
 
{
335
 
  void *boundminx=NULL,*boundmaxx=NULL,*wolfeminx=NULL;
336
 
  double boundmind=0,boundmaxd=0,wolfemind=0,aboundmind=0,aboundmaxd=0,awolfemind=0;
337
 
  long boundminl=0,boundmaxl=0,wolfeminl=0,aboundminl=0,aboundmaxl=0,awolfeminl=0;
338
 
  float boundminf=0,boundmaxf=0,wolfeminf=0,aboundminf=0,aboundmaxf=0,awolfeminf=0;
339
 
  int boundmini=0,boundmaxi=0,wolfemini=0,aboundmini=0,aboundmaxi=0,awolfemini=0;
340
 
 
341
 
  int g_a, g_b, g_c, g_d, g_e;
342
 
  int g_f, g_g, g_h, g_i, g_j;
343
 
  int g_k, g_l, g_m, g_n;
344
 
  int n = N;
345
 
  int me = GA_Nodeid ();
346
 
  int i;
347
 
  int dims[MAXDIM];
348
 
  int lo[MAXDIM], hi[MAXDIM];
349
 
  int index[MAXDIM];
350
 
  int index2[MAXDIM];
351
 
  int index3[MAXDIM];
352
 
  int block_size[MAXDIM], proc_grid[MAXDIM], proc_cnt;
353
 
  int needs_scaled_result;
354
 
  void *val=NULL;
355
 
  void *val3=NULL;
356
 
  int ival = -2;
357
 
  double dval = -2.0;
358
 
  float fval = -2.0;
359
 
  long lval = -2;
360
 
  DoubleComplex dcval;
361
 
  SingleComplex fcval;
362
 
  void *val2=NULL;
363
 
  void *val4=NULL;
364
 
  int ival2 = -3;
365
 
  double dval2 = -3.0;
366
 
  float fval2 = -3.0;
367
 
  long lval2 = -3;
368
 
  void *val5=NULL;
369
 
  int ival5 = 5;
370
 
  long lval5 = 5;
371
 
  double dval5 = 5.0;
372
 
  float fval5 = 5.0;
373
 
  DoubleComplex dcval2;
374
 
  DoubleComplex dcval3;
375
 
  DoubleComplex dcval4;
376
 
  DoubleComplex dcval5;
377
 
  DoubleComplex dcval6;
378
 
  DoubleComplex dcval7;
379
 
  SingleComplex fcval2;
380
 
  SingleComplex fcval3;
381
 
  SingleComplex fcval4;
382
 
  SingleComplex fcval5;
383
 
  SingleComplex fcval6;
384
 
  SingleComplex fcval7;
385
 
  void *vresult=NULL;
386
 
  int ivresult;
387
 
  double dvresult;
388
 
  float fvresult;
389
 
  long lvresult;
390
 
  DoubleComplex dcvresult;
391
 
  SingleComplex fcvresult;
392
 
  void *bvresult=NULL;
393
 
  DoubleComplex dcbvresult;
394
 
  SingleComplex fcbvresult;
395
 
  int ok = 1;
396
 
  int result=0,result2=0,result3=0;
397
 
  void *max=NULL;
398
 
  int imax;
399
 
  float fmax;
400
 
  long lmax;
401
 
  double dmax;
402
 
  DoubleComplex dcmax;
403
 
  SingleComplex fcmax;
404
 
  void *max2=NULL;
405
 
  DoubleComplex dcmax2;
406
 
  SingleComplex fcmax2;
407
 
  void *max3=NULL;
408
 
  DoubleComplex dcmax3;
409
 
  SingleComplex fcmax3;
410
 
 
411
 
  void *alpha=NULL, *beta=NULL;
412
 
  int ai = 1, bi = -1;
413
 
  long al = 1, bl = -1;
414
 
  float af = 1.0, bf = -1.0;
415
 
  double ad = 1.0, bd = -1.0;
416
 
  DoubleComplex adc, bdc;
417
 
  SingleComplex afc, bfc;
418
 
  double x1, x2, x3, x4;
419
 
  float fx1, fx2, fx3, fx4;
420
 
  void    *resultx=NULL;
421
 
  long    resultl=0,aresultl=0;
422
 
  double  resultd=0,aresultd=0;
423
 
  float   resultf=0,aresultf=0;
424
 
  int resulti=0,aresulti=0;
425
 
 
426
 
  adc.real = 1.0;
427
 
  adc.imag = 0.0;
428
 
  bdc.real = -1.0;
429
 
  bdc.imag = 0.0;
430
 
 
431
 
  afc.real = 1.0;
432
 
  afc.imag = 0.0;
433
 
  bfc.real = -1.0;
434
 
  bfc.imag = 0.0;
435
 
 
436
 
  needs_scaled_result = 0;
437
 
  
438
 
  dcval.real = -sin (3.0);
439
 
  dcval.imag = -cos (3.0);
440
 
  dcval2.real = 2 * sin (3.0);
441
 
  dcval2.imag = 2 * cos (3.0);
442
 
  dcval3.real = dcval.real*1.0e200;
443
 
  dcval3.imag = dcval.imag*1.0e200;
444
 
  dcval4.real = dcval2.real*1.0e200;
445
 
  dcval4.imag = dcval2.imag*1.0e200;
446
 
  dcval5.real = 5.0;
447
 
  dcval5.imag = 0.0;
448
 
  dcval6.real = dcval3.imag;
449
 
  dcval6.imag = dcval3.real;
450
 
  dcval7.real = dcval4.imag;
451
 
  dcval7.imag = dcval4.real;
452
 
  
453
 
  fcval.real = -sin (3.0);
454
 
  fcval.imag = -cos (3.0);
455
 
  fcval2.real = 2 * sin (3.0);
456
 
  fcval2.imag = 2 * cos (3.0);
457
 
  fcval3.real = fcval.real*1.0e200;
458
 
  fcval3.imag = fcval.imag*1.0e200;
459
 
  fcval4.real = fcval2.real*1.0e200;
460
 
  fcval4.imag = fcval2.imag*1.0e200;
461
 
  fcval5.real = 5.0;
462
 
  fcval5.imag = 0.0;
463
 
  fcval6.real = fcval3.imag;
464
 
  fcval6.imag = fcval3.real;
465
 
  fcval7.real = fcval4.imag;
466
 
  fcval7.imag = fcval4.real;
467
 
  
468
 
  
469
 
  proc_cnt=1;
470
 
  for (i = 0; i < dim; i++) {
471
 
     dims[i] = N;
472
 
     block_size[i] = BLOCK_SIZE;
473
 
     if (i<dim-1 && proc_cnt < GA_Nnodes()) {
474
 
        proc_grid[i] = 2;
475
 
        proc_cnt *= 2;
476
 
     } else if (proc_cnt >= GA_Nnodes()) {
477
 
        proc_grid[i] = 1;
478
 
     } else {
479
 
        proc_grid[i] = GA_Nnodes()/proc_cnt;
480
 
     }
481
 
  }
482
 
 
483
 
  for (i = 0; i < dim; i++)
484
 
  {
485
 
     lo[i] = 0;
486
 
     hi[i] = N - 1;
487
 
  }
488
 
#if BLOCK_CYCLIC
489
 
  g_a = GA_Create_handle();
490
 
  GA_Set_data(g_a,dim,dims,type);
491
 
  GA_Set_array_name(g_a,"A");
492
 
# if USE_SCALAPACK
493
 
  GA_Set_block_cyclic_proc_grid(g_a,block_size,proc_grid);
494
 
# else
495
 
  GA_Set_block_cyclic(g_a,block_size);
496
 
# endif
497
 
  GA_Allocate(g_a);
498
 
#else
499
 
  g_a = NGA_Create (type, dim, dims, "A", NULL);
500
 
  if (!g_a)
501
 
    GA_Error ("create failed: A", n);
502
 
#endif
503
 
 
504
 
  g_b = GA_Duplicate (g_a, "B");
505
 
  if (!g_b)
506
 
    GA_Error ("duplicate failed: B", n);
507
 
 
508
 
  g_c = GA_Duplicate (g_a, "C");
509
 
  if (!g_c)
510
 
    GA_Error ("duplicate failed: C", n);
511
 
 
512
 
  g_d = GA_Duplicate (g_a, "D");
513
 
  if (!g_d)
514
 
    GA_Error ("duplicate failed: D", n);
515
 
 
516
 
  g_e = GA_Duplicate (g_a, "E");
517
 
  if (!g_e)
518
 
    GA_Error ("duplicate failed: E", n);
519
 
 
520
 
  g_f = GA_Duplicate (g_a, "F");
521
 
  if (!g_f)
522
 
    GA_Error ("duplicate failed: F", n);
523
 
 
524
 
  g_g = GA_Duplicate (g_a, "G");
525
 
  if (!g_g)
526
 
    GA_Error ("duplicate failed: G", n);
527
 
 
528
 
  g_h = GA_Duplicate (g_a, "H");
529
 
  if (!g_h)
530
 
    GA_Error ("duplicate failed: H", n);
531
 
 
532
 
  g_i = GA_Duplicate (g_a, "I");
533
 
  if (!g_i)
534
 
    GA_Error ("duplicate failed: I", n);
535
 
 
536
 
  g_j = GA_Duplicate (g_a, "J");
537
 
  if (!g_j)
538
 
    GA_Error ("duplicate failed: J", n);
539
 
 
540
 
  g_k = GA_Duplicate (g_a, "K");
541
 
  if (!g_k)
542
 
    GA_Error ("duplicate failed: K", n);
543
 
 
544
 
  g_l = GA_Duplicate (g_a, "L");
545
 
  if (!g_l)
546
 
    GA_Error ("duplicate failed: L", n);
547
 
 
548
 
  g_m = GA_Duplicate (g_a, "M");
549
 
  if (!g_m)
550
 
    GA_Error ("duplicate failed: M", n);
551
 
 
552
 
  g_n = GA_Duplicate (g_a, "N");
553
 
  if (!g_m)
554
 
    GA_Error ("duplicate failed: N", n);
555
 
  /*initialize  with zero */
556
 
  GA_Zero (g_a);
557
 
  GA_Zero (g_b);
558
 
  GA_Zero (g_c);
559
 
  GA_Zero (g_d);
560
 
  GA_Zero (g_e);
561
 
  GA_Zero (g_f);
562
 
  GA_Zero (g_g);
563
 
  GA_Zero (g_h);
564
 
  GA_Zero (g_i);
565
 
  GA_Zero (g_j);
566
 
  GA_Zero (g_k);
567
 
  GA_Zero (g_l);
568
 
  GA_Zero (g_m);
569
 
  GA_Zero (g_n);
570
 
 
571
 
  switch (type)
572
 
    {
573
 
    case C_INT:
574
 
      val = &ival;
575
 
      val2 = &ival2;
576
 
      val5 = &ival5;
577
 
      vresult = &ivresult;
578
 
      resultx = &resulti;
579
 
      boundminx = &boundmini;
580
 
      boundmaxx = &boundmaxi;
581
 
      wolfeminx = &wolfemini;
582
 
      break;
583
 
    case C_DCPL:
584
 
      val = &dcval;
585
 
      val2 = &dcval2;
586
 
      val3 = &dcval3;
587
 
      val4 = &dcval4;
588
 
      val5 = &dcval5;
589
 
      vresult = &dcvresult;
590
 
      bvresult = &dcbvresult;
591
 
      break;
592
 
    case C_SCPL:
593
 
      val = &fcval;
594
 
      val2 = &fcval2;
595
 
      val3 = &fcval3;
596
 
      val4 = &fcval4;
597
 
      val5 = &fcval5;
598
 
      vresult = &fcvresult;
599
 
      bvresult = &fcbvresult;
600
 
      break;
601
 
 
602
 
    case C_DBL:
603
 
      val = &dval;
604
 
      val2 = &dval2;
605
 
      val5 = &dval5;
606
 
      vresult = &dvresult;
607
 
      resultx = &resultd;
608
 
      boundminx = &boundmind;
609
 
      boundmaxx = &boundmaxd;
610
 
      wolfeminx = &wolfemind;
611
 
      break;
612
 
    case C_FLOAT:
613
 
      val = &fval;
614
 
      val2 = &fval2;
615
 
      val5 = &fval5;
616
 
      vresult = &fvresult;
617
 
      resultx = &resultf;
618
 
      boundminx = &boundminf;
619
 
      boundmaxx = &boundmaxf;
620
 
      wolfeminx = &wolfeminf;
621
 
      break;
622
 
    case C_LONG:
623
 
      val = &lval;
624
 
      val2 = &lval2;
625
 
      val5 = &lval5;
626
 
      vresult = &lvresult;
627
 
      resultx = &resultl;
628
 
      boundminx = &boundminl;
629
 
      boundmaxx = &boundmaxl;
630
 
      wolfeminx = &wolfeminl;
631
 
      break;
632
 
    default:
633
 
      GA_Error ("wrong data type.", type);
634
 
    }
635
 
 
636
 
 
637
 
  NGA_Fill_patch (g_a, lo, hi, val);
638
 
  NGA_Fill_patch (g_b, lo, hi, val2);
639
 
  NGA_Pnfill_patch (g_j, lo, hi);
640
 
  switch (OP)
641
 
    {
642
 
      double tmp, tmp2;
643
 
      float  tmpf, tmp2f;
644
 
    case OP_ABS:
645
 
      if (me == 0)
646
 
    printf ("Testing GA_Abs_value...");
647
 
      GA_Abs_value_patch (g_a, lo, hi);
648
 
      ivresult = GA_ABS (ival);
649
 
      dvresult = GA_ABS (dval);
650
 
      fvresult = GA_ABS (fval);
651
 
      lvresult = GA_ABS (lval);
652
 
      
653
 
      if (GA_ABS(dcval.real) >= GA_ABS(dcval.imag)) {
654
 
    if (dcval.real == (double)0.0) {
655
 
      dcvresult.real = (double)0.0;
656
 
    } else {
657
 
      x1 = dcval.imag/dcval.real;
658
 
      dcvresult.real = GA_ABS(dcval.real)*sqrt(((double)1.0)+(x1*x1));
659
 
    }
660
 
      } else {
661
 
    x1 = dcval.real/dcval.imag;
662
 
    dcvresult.real = GA_ABS(dcval.imag)*sqrt(((double)1.0)+(x1*x1));
663
 
      }
664
 
      dcvresult.imag = 0.0;
665
 
      NGA_Fill_patch (g_d, lo, hi, vresult);
666
 
      
667
 
      if (type == C_DCPL) {
668
 
        needs_scaled_result = 1;
669
 
    NGA_Fill_patch(g_f,lo,hi,val3);
670
 
    GA_Abs_value_patch (g_f, lo, hi);
671
 
    if (GA_ABS(dcval3.real) >= GA_ABS(dcval3.imag)) {
672
 
      if (dcval3.real == (double)0.0) {
673
 
        dcbvresult.real = (double)0.0;
674
 
      } else {
675
 
        x1 = dcval3.imag/dcval3.real;
676
 
        dcbvresult.real = GA_ABS(dcval3.real)*sqrt(((double)1.0)+(x1*x1));
677
 
      }
678
 
    } else {
679
 
      x1 = dcval3.real/dcval3.imag;
680
 
      dcbvresult.real = GA_ABS(dcval3.imag)*sqrt(((double)1.0)+(x1*x1));
681
 
    }
682
 
    dcbvresult.imag = (double)0.0;
683
 
    NGA_Fill_patch (g_i, lo, hi, bvresult);
684
 
    NGA_Fill_patch(g_k,lo,hi,&dcval6);
685
 
    GA_Abs_value_patch (g_k, lo, hi);
686
 
    if (GA_ABS(dcval6.real) >= GA_ABS(dcval6.imag)) {
687
 
      if (dcval6.real == (double)0.0) {
688
 
        dcbvresult.real = (double)0.0;
689
 
      } else {
690
 
        x1 = dcval6.imag/dcval6.real;
691
 
        dcbvresult.real = GA_ABS(dcval6.real)*sqrt(((double)1.0)+(x1*x1));
692
 
      }
693
 
    } else {
694
 
      x1 = dcval6.real/dcval6.imag;
695
 
      dcbvresult.real = GA_ABS(dcval6.imag)*sqrt(((double)1.0)+(x1*x1));
696
 
    }
697
 
    NGA_Fill_patch (g_n, lo, hi, bvresult);
698
 
      }
699
 
      if (type == C_SCPL) {
700
 
        needs_scaled_result = 1;
701
 
    NGA_Fill_patch(g_f,lo,hi,val3);
702
 
    GA_Abs_value_patch (g_f, lo, hi);
703
 
    if (GA_ABS(fcval3.real) >= GA_ABS(fcval3.imag)) {
704
 
      if (fcval3.real == (float )0.0) {
705
 
        fcbvresult.real = (float )0.0;
706
 
      } else {
707
 
        fx1 = fcval3.imag/fcval3.real;
708
 
        fcbvresult.real = GA_ABS(fcval3.real)*sqrt(((float)1.0)+(fx1*fx1));
709
 
      }
710
 
    } else {
711
 
      fx1 = fcval3.real/fcval3.imag;
712
 
      fcbvresult.real = GA_ABS(fcval3.imag)*sqrt(((float)1.0)+(fx1*fx1));
713
 
    }
714
 
    fcbvresult.imag = (float)0.0;
715
 
    NGA_Fill_patch (g_i, lo, hi, bvresult);
716
 
    NGA_Fill_patch(g_k,lo,hi,&dcval6);
717
 
    GA_Abs_value_patch (g_k, lo, hi);
718
 
    if (GA_ABS(fcval6.real) >= GA_ABS(fcval6.imag)) {
719
 
      if (fcval6.real == (float)0.0) {
720
 
        fcbvresult.real = (float)0.0;
721
 
      } else {
722
 
        fx1 = fcval6.imag/fcval6.real;
723
 
        fcbvresult.real = GA_ABS(fcval6.real)*sqrt(((float)1.0)+(fx1*fx1));
724
 
      }
725
 
    } else {
726
 
      fx1 = fcval6.real/fcval6.imag;
727
 
      fcbvresult.real = GA_ABS(fcval6.imag)*sqrt(((float)1.0)+(fx1*fx1));
728
 
    }
729
 
    NGA_Fill_patch (g_n, lo, hi, bvresult);
730
 
      }
731
 
      break;
732
 
    case OP_ADD_CONST:
733
 
      if (me == 0)
734
 
    printf ("Testing GA_Add_const...");
735
 
      GA_Add_constant_patch (g_a, lo, hi, val2);
736
 
      ivresult = ival + ival2;
737
 
      dvresult = dval + dval2;
738
 
      fvresult = fval + fval2;
739
 
      lvresult = lval + lval2;
740
 
      dcvresult.real = dcval.real + dcval2.real;
741
 
      dcvresult.imag = dcval.imag + dcval2.imag;
742
 
      fcvresult.real = fcval.real + fcval2.real;
743
 
      fcvresult.imag = fcval.imag + fcval2.imag;
744
 
      NGA_Fill_patch (g_d, lo, hi, vresult);
745
 
      break;
746
 
    case OP_RECIP:
747
 
      if (me == 0)
748
 
    printf ("Testing GA_Recip...");
749
 
      GA_Recip_patch (g_a, lo, hi);
750
 
      ivresult = ((int)1) / ival;
751
 
      dvresult = ((double)1.0) / dval;
752
 
      fvresult = ((float)1.0) / fval;
753
 
      lvresult = ((long)1) / lval;
754
 
      
755
 
      if (GA_ABS(dcval.real) >= GA_ABS(dcval.imag)) {
756
 
    if (dcval.real != (double)0.0) {
757
 
      tmp = dcval.imag/dcval.real;
758
 
      tmp2 = ((double)1.0)/((((double)1.0)+(tmp*tmp))*dcval.real);
759
 
      dcvresult.real = tmp2;
760
 
      dcvresult.imag = -tmp * tmp2;
761
 
    } else {
762
 
      printf("Error in testing GA_Recip dcval = 0.0\n");
763
 
    }
764
 
      } else {
765
 
    tmp = dcval.real/dcval.imag;
766
 
    tmp2 = ((double)1.0)/((((double)1.0)+(tmp*tmp))*dcval.imag);
767
 
    dcvresult.real = tmp * tmp2;
768
 
    dcvresult.imag = -tmp2;
769
 
      }
770
 
      NGA_Fill_patch (g_d, lo, hi, vresult);
771
 
      
772
 
      if (type == C_DCPL) {
773
 
        needs_scaled_result = 1;
774
 
    NGA_Fill_patch (g_f, lo, hi, val3);
775
 
    GA_Recip_patch (g_f, lo, hi);
776
 
    if (GA_ABS(dcval3.real) >= GA_ABS(dcval3.imag)) {
777
 
      if (dcval3.real == (double)0.0) {
778
 
        printf("Error testing GA_Recip, dcval3.real = 0.0\n");
779
 
      } else {
780
 
        tmp = dcval3.imag/dcval3.real;
781
 
        tmp2 = ((double)1.0)/((((double)1.0)+(tmp*tmp))*dcval3.real);
782
 
            dcbvresult.real = tmp2;
783
 
        dcbvresult.imag = -tmp * tmp2;
784
 
      }
785
 
    } else {
786
 
      tmp = dcval3.real/dcval3.imag;
787
 
      tmp2 = ((double)1.0)/((((double)1.0)+(tmp*tmp))*dcval3.imag);
788
 
      dcbvresult.real = tmp * tmp2;
789
 
      dcbvresult.imag = -tmp2;
790
 
    }
791
 
    NGA_Fill_patch (g_i, lo, hi, bvresult);
792
 
    NGA_Fill_patch(g_k,lo,hi,&dcval6);
793
 
    GA_Recip_patch (g_k, lo, hi);
794
 
    if (GA_ABS(dcval6.real) >= GA_ABS(dcval6.imag)) {
795
 
      if (dcval6.real == (double)0.0) {
796
 
        printf("Error testing GA_Recip, dcval6.real = 0.0\n");
797
 
      } else {
798
 
        tmp = dcval6.imag/dcval6.real;
799
 
        tmp2 = ((double)1.0)/((((double)1.0)+(tmp*tmp))*dcval6.real);
800
 
            dcbvresult.real = tmp2;
801
 
        dcbvresult.imag = -tmp * tmp2;
802
 
      }
803
 
    } else {
804
 
      tmp = dcval6.real/dcval6.imag;
805
 
      tmp2 = ((double)1.0)/((((double)1.0)+(tmp*tmp))*dcval6.imag);
806
 
      dcbvresult.real = tmp * tmp2;
807
 
      dcbvresult.imag = -tmp2;
808
 
    }
809
 
    NGA_Fill_patch (g_n, lo, hi, bvresult);
810
 
      }
811
 
      if (type == C_SCPL) {
812
 
        needs_scaled_result = 1;
813
 
    NGA_Fill_patch (g_f, lo, hi, val3);
814
 
    GA_Recip_patch (g_f, lo, hi);
815
 
    if (GA_ABS(fcval3.real) >= GA_ABS(fcval3.imag)) {
816
 
      if (fcval3.real == (float )0.0) {
817
 
        printf("Error testing GA_Recip, fcval3.real = 0.0\n");
818
 
      } else {
819
 
        tmpf = fcval3.imag/fcval3.real;
820
 
        tmp2f = ((float)1.0)/((((float)1.0)+(tmpf*tmpf))*fcval3.real);
821
 
            fcbvresult.real = tmp2f;
822
 
        fcbvresult.imag = -tmpf * tmp2f;
823
 
      }
824
 
    } else {
825
 
      tmpf = fcval3.real/fcval3.imag;
826
 
      tmp2f = ((float)1.0)/((((float)1.0)+(tmpf*tmpf))*fcval3.imag);
827
 
      fcbvresult.real = tmpf * tmp2f;
828
 
      fcbvresult.imag = -tmp2f;
829
 
    }
830
 
    NGA_Fill_patch (g_i, lo, hi, bvresult);
831
 
    NGA_Fill_patch(g_k,lo,hi,&dcval6);
832
 
    GA_Recip_patch (g_k, lo, hi);
833
 
    if (GA_ABS(fcval6.real) >= GA_ABS(fcval6.imag)) {
834
 
      if (fcval6.real == (float)0.0) {
835
 
        printf("Error testing GA_Recip, fcval6.real = 0.0\n");
836
 
      } else {
837
 
        tmpf = fcval6.imag/fcval6.real;
838
 
        tmp2f = ((float)1.0)/((((float)1.0)+(tmpf*tmpf))*fcval6.real);
839
 
            fcbvresult.real = tmp2f;
840
 
        fcbvresult.imag = -tmpf * tmp2f;
841
 
      }
842
 
    } else {
843
 
      tmpf = fcval6.real/fcval6.imag;
844
 
      tmp2f = ((float)1.0)/((((float)1.0)+(tmpf*tmpf))*fcval6.imag);
845
 
      fcbvresult.real = tmpf * tmp2f;
846
 
      fcbvresult.imag = -tmp2f;
847
 
    }
848
 
    NGA_Fill_patch (g_n, lo, hi, bvresult);
849
 
      }
850
 
      break;
851
 
    case OP_ELEM_MULT:
852
 
      if (me == 0)
853
 
    printf ("Testing GA_Elem_multiply...");
854
 
      NGA_Fill_patch (g_b, lo, hi, val2);
855
 
      /* g_c is different from g_a or g_b*/
856
 
      GA_Elem_multiply_patch (g_a, lo, hi, g_b, lo, hi, g_c, lo, hi);
857
 
 
858
 
      ivresult = ival * ival2;
859
 
      dvresult = dval * dval2;
860
 
      fvresult = fval * fval2;
861
 
      lvresult = lval * lval2;
862
 
      dcvresult.real = dcval.real * dcval2.real - dcval.imag * dcval2.imag;
863
 
      dcvresult.imag = dcval.real * dcval2.imag + dcval2.real * dcval.imag;
864
 
      NGA_Fill_patch (g_d, lo, hi, vresult);
865
 
      break;
866
 
    case OP_ELEM_DIV:
867
 
      if (me == 0)
868
 
    printf ("Testing GA_Elem_divide...");
869
 
      NGA_Fill_patch (g_b, lo, hi, val2);
870
 
      GA_Elem_divide_patch (g_a, lo, hi, g_b, lo, hi, g_c, lo, hi);
871
 
      ivresult = ival / ival2;
872
 
      dvresult = dval / dval2;
873
 
      fvresult = fval / fval2;
874
 
      lvresult = lval / lval2;
875
 
      dcvresult.real = 0.0;
876
 
      dcvresult.imag = 0.0;
877
 
      
878
 
      if (GA_ABS(dcval2.real) >= GA_ABS(dcval2.imag)) {
879
 
    if (dcval2.real != (double)0.0) {
880
 
      tmp = dcval2.imag/dcval2.real;
881
 
      tmp2 = ((double)1.0)/(dcval2.real*(((double)1.0)+(tmp*tmp)));
882
 
      dcvresult.real = (dcval.real + dcval.imag*tmp)*tmp2;
883
 
      dcvresult.imag = (dcval.imag - dcval.real*tmp)*tmp2;
884
 
    } else {
885
 
      printf("Error in testing GA_Elem_divide dcval = 0.0\n");
886
 
    } 
887
 
      } else {
888
 
    tmp = dcval2.real/dcval2.imag;
889
 
    tmp2 = 1.0/(dcval2.imag*(1.0+(tmp*tmp)));
890
 
    dcvresult.real = (dcval.real*tmp + dcval.imag)*tmp2;
891
 
    dcvresult.imag = (dcval.imag*tmp - dcval.real)*tmp2;
892
 
      }
893
 
      NGA_Fill_patch (g_d, lo, hi, vresult);
894
 
      
895
 
      if (type == C_DCPL) {
896
 
        needs_scaled_result = 1;
897
 
    NGA_Fill_patch (g_f, lo, hi, val3);
898
 
    NGA_Fill_patch (g_g, lo, hi, val4);
899
 
    GA_Elem_divide_patch (g_f, lo, hi, g_g, lo, hi, g_h, lo, hi);
900
 
    dcbvresult.real = (double)0.0;
901
 
    dcbvresult.imag = (double)0.0;
902
 
    if (GA_ABS(dcval4.real) >= GA_ABS(dcval4.imag)) {
903
 
      if (dcval4.real != (double)0.0) {
904
 
        tmp = dcval4.imag/dcval4.real;
905
 
        tmp2 = ((double)1.0)/(dcval4.real*(((double)1.0)+(tmp*tmp)));
906
 
        dcbvresult.real = (dcval3.real + dcval3.imag*tmp)*tmp2;
907
 
        dcbvresult.imag = (dcval3.imag - dcval3.real*tmp)*tmp2;
908
 
      } else {
909
 
        printf("Error in testing GA_Elem_divide dcval4 = 0.0\n");
910
 
      } 
911
 
    } else {
912
 
      tmp = dcval4.real/dcval4.imag;
913
 
      tmp2 = ((double)1.0)/(dcval4.imag*(((double)1.0)+(tmp*tmp)));
914
 
      dcbvresult.real = (dcval3.real*tmp + dcval3.imag)*tmp2;
915
 
      dcbvresult.imag = (dcval3.imag*tmp - dcval3.real)*tmp2;
916
 
    }
917
 
    NGA_Fill_patch (g_i, lo, hi, bvresult);
918
 
    NGA_Fill_patch (g_k, lo, hi, &dcval6);
919
 
    NGA_Fill_patch (g_l, lo, hi, &dcval7);
920
 
    GA_Elem_divide_patch (g_k, lo, hi, g_l, lo, hi, g_m, lo, hi);
921
 
    dcbvresult.real = (double)0.0;
922
 
    dcbvresult.imag = (double)0.0;
923
 
    if (GA_ABS(dcval7.real) >= GA_ABS(dcval7.imag)) {
924
 
      if (dcval7.real != (double)0.0) {
925
 
        tmp = dcval7.imag/dcval7.real;
926
 
        tmp2 = ((double)1.0)/(dcval7.real*(((double)1.0)+(tmp*tmp)));
927
 
        dcbvresult.real = (dcval6.real + dcval6.imag*tmp)*tmp2;
928
 
        dcbvresult.imag = (dcval6.imag - dcval6.real*tmp)*tmp2;
929
 
      } else {
930
 
        printf("Error in testing GA_Elem_divide dcval7 = 0.0\n");
931
 
      } 
932
 
    } else {
933
 
      tmp = dcval7.real/dcval7.imag;
934
 
      tmp2 = ((double)1.0)/(dcval7.imag*(((double)1.0)+(tmp*tmp)));
935
 
      dcbvresult.real = (dcval6.real*tmp + dcval6.imag)*tmp2;
936
 
      dcbvresult.imag = (dcval6.imag*tmp - dcval6.real)*tmp2;
937
 
    }
938
 
    NGA_Fill_patch (g_n, lo, hi, bvresult);
939
 
      }
940
 
      if (type == C_SCPL) {
941
 
        needs_scaled_result = 1;
942
 
    NGA_Fill_patch (g_f, lo, hi, val3);
943
 
    NGA_Fill_patch (g_g, lo, hi, val4);
944
 
    GA_Elem_divide_patch (g_f, lo, hi, g_g, lo, hi, g_h, lo, hi);
945
 
    fcbvresult.real = (float)0.0;
946
 
    fcbvresult.imag = (float)0.0;
947
 
    if (GA_ABS(fcval4.real) >= GA_ABS(fcval4.imag)) {
948
 
      if (fcval4.real != (float)0.0) {
949
 
        tmpf = fcval4.imag/fcval4.real;
950
 
        tmp2f = ((float)1.0)/(fcval4.real*(((float)1.0)+(tmpf*tmpf)));
951
 
        fcbvresult.real = (fcval3.real + fcval3.imag*tmpf)*tmp2f;
952
 
        fcbvresult.imag = (fcval3.imag - fcval3.real*tmpf)*tmp2f;
953
 
      } else {
954
 
        printf("Error in testing GA_Elem_divide fcval4 = 0.0\n");
955
 
      } 
956
 
    } else {
957
 
      tmpf = fcval4.real/fcval4.imag;
958
 
      tmp2f = ((float)1.0)/(fcval4.imag*(((float)1.0)+(tmpf*tmpf)));
959
 
      fcbvresult.real = (fcval3.real*tmpf + fcval3.imag)*tmp2f;
960
 
      fcbvresult.imag = (fcval3.imag*tmpf - fcval3.real)*tmp2f;
961
 
    }
962
 
    NGA_Fill_patch (g_i, lo, hi, bvresult);
963
 
    NGA_Fill_patch (g_k, lo, hi, &dcval6);
964
 
    NGA_Fill_patch (g_l, lo, hi, &dcval7);
965
 
    GA_Elem_divide_patch (g_k, lo, hi, g_l, lo, hi, g_m, lo, hi);
966
 
    fcbvresult.real = (float)0.0;
967
 
    fcbvresult.imag = (float)0.0;
968
 
    if (GA_ABS(fcval7.real) >= GA_ABS(fcval7.imag)) {
969
 
      if (fcval7.real != (float)0.0) {
970
 
        tmpf = fcval7.imag/fcval7.real;
971
 
        tmp2f = ((float)1.0)/(fcval7.real*(((float)1.0)+(tmpf*tmpf)));
972
 
        fcbvresult.real = (fcval6.real + fcval6.imag*tmpf)*tmp2f;
973
 
        fcbvresult.imag = (fcval6.imag - fcval6.real*tmpf)*tmp2f;
974
 
      } else {
975
 
        printf("Error in testing GA_Elem_divide fcval7 = 0.0\n");
976
 
      } 
977
 
    } else {
978
 
      tmpf = fcval7.real/fcval7.imag;
979
 
      tmp2f = ((float)1.0)/(fcval7.imag*(((float)1.0)+(tmpf*tmpf)));
980
 
      fcbvresult.real = (fcval6.real*tmpf + fcval6.imag)*tmp2f;
981
 
      fcbvresult.imag = (fcval6.imag*tmpf - fcval6.real)*tmp2f;
982
 
    }
983
 
    NGA_Fill_patch (g_n, lo, hi, bvresult);
984
 
      }
985
 
      break;
986
 
 
987
 
    case OP_ELEM_MAX:
988
 
      if (me == 0)
989
 
    printf ("Testing GA_Elem_maximum...");
990
 
      /*NGA_Fill_patch (g_b, lo, hi, val2);*/
991
 
      GA_Elem_maximum_patch (g_a, lo, hi, g_b, lo, hi, g_c, lo, hi);
992
 
      ivresult = GA_MAX (ival, ival2);
993
 
      dvresult = GA_MAX (dval, dval2);
994
 
      fvresult = GA_MAX (fval, fval2);
995
 
      lvresult = GA_MAX (lval, lval2);
996
 
      tmp  = GA_MAX(GA_ABS(dcval.real),GA_ABS(dcval.imag));
997
 
      tmp2 = GA_MAX(GA_ABS(dcval2.real),GA_ABS(dcval2.imag));
998
 
      tmp  = GA_MAX(tmp,tmp2);
999
 
      dcvresult.real = dcval.real;
1000
 
      dcvresult.imag = dcval.imag;
1001
 
      if (tmp != 0.0) {
1002
 
    tmp = ((double)1.0)/tmp;
1003
 
    x1 = dcval.real*tmp;
1004
 
    x2 = dcval.imag*tmp;
1005
 
    x3 = dcval2.real*tmp;
1006
 
    x4 = dcval2.imag*tmp;
1007
 
    tmp = x1*x1 + x2*x2;
1008
 
    tmp2 = x3*x3 + x4*x4;
1009
 
    if (tmp2 > tmp) {
1010
 
      dcvresult.real = dcval2.real;
1011
 
      dcvresult.imag = dcval2.imag;
1012
 
    }
1013
 
      }
1014
 
      NGA_Fill_patch (g_d, lo, hi, vresult);
1015
 
      if (type == C_DCPL) {
1016
 
        needs_scaled_result = 1;
1017
 
    NGA_Fill_patch (g_f, lo, hi, val3);
1018
 
    NGA_Fill_patch (g_g, lo, hi, val4);
1019
 
    GA_Elem_maximum_patch (g_f, lo, hi, g_g, lo, hi, g_h, lo, hi);
1020
 
    tmp  = GA_MAX(GA_ABS(dcval3.real),GA_ABS(dcval3.imag));
1021
 
    tmp2 = GA_MAX(GA_ABS(dcval4.real),GA_ABS(dcval4.imag));
1022
 
    tmp  = GA_MAX(tmp,tmp2);
1023
 
    dcvresult.real = dcval3.real;
1024
 
    dcvresult.imag = dcval3.imag;
1025
 
    if (tmp != 0.0) {
1026
 
      tmp = ((double)1.0)/tmp;
1027
 
      x1 = dcval3.real*tmp;
1028
 
      x2 = dcval3.imag*tmp;
1029
 
      x3 = dcval4.real*tmp;
1030
 
      x4 = dcval4.imag*tmp;
1031
 
      tmp = x1*x1 + x2*x2;
1032
 
      tmp2 = x3*x3 + x4*x4;
1033
 
      if (tmp2 > tmp) {
1034
 
        dcvresult.real = dcval4.real;
1035
 
        dcvresult.imag = dcval4.imag;
1036
 
      }
1037
 
    }
1038
 
    NGA_Fill_patch (g_i, lo, hi, vresult);
1039
 
    NGA_Fill_patch (g_k, lo, hi, &dcval6);
1040
 
    NGA_Fill_patch (g_l, lo, hi, &dcval7);
1041
 
    GA_Elem_maximum_patch (g_k, lo, hi, g_l, lo, hi, g_m, lo, hi);
1042
 
    tmp  = GA_MAX(GA_ABS(dcval6.real),GA_ABS(dcval6.imag));
1043
 
    tmp2 = GA_MAX(GA_ABS(dcval7.real),GA_ABS(dcval7.imag));
1044
 
    tmp  = GA_MAX(tmp,tmp2);
1045
 
    dcvresult.real = dcval6.real;
1046
 
    dcvresult.imag = dcval6.imag;
1047
 
    if (tmp != 0.0) {
1048
 
      tmp = ((double)1.0)/tmp;
1049
 
      x1 = dcval6.real*tmp;
1050
 
      x2 = dcval6.imag*tmp;
1051
 
      x3 = dcval7.real*tmp;
1052
 
      x4 = dcval7.imag*tmp;
1053
 
      tmp = x1*x1 + x2*x2;
1054
 
      tmp2 = x3*x3 + x4*x4;
1055
 
      if (tmp2 > tmp) {
1056
 
        dcvresult.real = dcval7.real;
1057
 
        dcvresult.imag = dcval7.imag;
1058
 
      }
1059
 
    }
1060
 
    NGA_Fill_patch (g_n, lo, hi, vresult);
1061
 
      }
1062
 
      break;
1063
 
    case OP_ELEM_MIN:
1064
 
      if (me == 0)
1065
 
    printf ("Testing GA_Elem_minimum...");
1066
 
      NGA_Fill_patch (g_b, lo, hi, val2);
1067
 
      GA_Elem_minimum_patch (g_a, lo, hi, g_b, lo, hi, g_c, lo, hi);
1068
 
      ivresult = GA_MIN (ival, ival2);
1069
 
      dvresult = GA_MIN (dval, dval2);
1070
 
      fvresult = GA_MIN (fval, fval2);
1071
 
      lvresult = GA_MIN (lval, lval2);
1072
 
      tmp  = GA_MAX(GA_ABS(dcval.real),GA_ABS(dcval.imag));
1073
 
      tmp2 = GA_MAX(GA_ABS(dcval2.real),GA_ABS(dcval2.imag));
1074
 
      tmp  = GA_MAX(tmp,tmp2);
1075
 
      dcvresult.real = dcval.real;
1076
 
      dcvresult.imag = dcval.imag;
1077
 
      if (tmp != 0.0) {
1078
 
    tmp = 1.0/tmp;
1079
 
    x1 = dcval.real*tmp;
1080
 
    x2 = dcval.imag*tmp;
1081
 
    x3 = dcval2.real*tmp;
1082
 
    x4 = dcval2.imag*tmp;
1083
 
    tmp = x1*x1 + x2*x2;
1084
 
    tmp2 = x3*x3 + x4*x4;
1085
 
    if (tmp2 < tmp) {
1086
 
      dcvresult.real = dcval2.real;
1087
 
      dcvresult.imag = dcval2.imag;
1088
 
    }
1089
 
      }
1090
 
      NGA_Fill_patch (g_d, lo, hi, vresult);
1091
 
      if (type == C_SCPL) {
1092
 
        needs_scaled_result = 1;
1093
 
    NGA_Fill_patch (g_f, lo, hi, val3);
1094
 
    NGA_Fill_patch (g_g, lo, hi, val4);
1095
 
    GA_Elem_minimum_patch (g_f, lo, hi, g_g, lo, hi, g_h, lo, hi);
1096
 
    tmpf  = GA_MAX(GA_ABS(fcval3.real),GA_ABS(fcval3.imag));
1097
 
    tmp2f = GA_MAX(GA_ABS(fcval4.real),GA_ABS(fcval4.imag));
1098
 
    tmpf  = GA_MAX(tmpf,tmp2f);
1099
 
    fcvresult.real = fcval3.real;
1100
 
    fcvresult.imag = fcval3.imag;
1101
 
    if (tmpf != 0.0) {
1102
 
      tmpf = ((float)1.0)/tmpf;
1103
 
      fx1 = fcval3.real*tmpf;
1104
 
      fx2 = fcval3.imag*tmpf;
1105
 
      fx3 = fcval4.real*tmpf;
1106
 
      fx4 = fcval4.imag*tmpf;
1107
 
      tmpf = fx1*fx1 + fx2*fx2;
1108
 
      tmp2f = fx3*fx3 + fx4*fx4;
1109
 
      if (tmp2f < tmpf) {
1110
 
        fcvresult.real = fcval4.real;
1111
 
        fcvresult.imag = fcval4.imag;
1112
 
      }
1113
 
    }
1114
 
    NGA_Fill_patch (g_i, lo, hi, vresult);
1115
 
    NGA_Fill_patch (g_k, lo, hi, &dcval6);
1116
 
    NGA_Fill_patch (g_l, lo, hi, &dcval7);
1117
 
    GA_Elem_minimum_patch (g_k, lo, hi, g_l, lo, hi, g_m, lo, hi);
1118
 
    tmpf  = GA_MAX(GA_ABS(fcval6.real),GA_ABS(fcval6.imag));
1119
 
    tmp2f = GA_MAX(GA_ABS(fcval7.real),GA_ABS(fcval7.imag));
1120
 
    tmpf  = GA_MAX(tmpf,tmp2f);
1121
 
    fcvresult.real = fcval6.real;
1122
 
    fcvresult.imag = fcval6.imag;
1123
 
    if (tmpf != 0.0) {
1124
 
      tmpf = ((float)1.0)/tmpf;
1125
 
      fx1 = fcval6.real*tmpf;
1126
 
      fx2 = fcval6.imag*tmpf;
1127
 
      fx3 = fcval7.real*tmpf;
1128
 
      fx4 = fcval7.imag*tmpf;
1129
 
      tmpf = fx1*fx1 + fx2*fx2;
1130
 
      tmp2f = fx3*fx3 + fx4*fx4;
1131
 
      if (tmp2f < tmpf) {
1132
 
        fcvresult.real = fcval7.real;
1133
 
        fcvresult.imag = fcval7.imag;
1134
 
      }
1135
 
    }
1136
 
    NGA_Fill_patch (g_n, lo, hi, vresult);
1137
 
      }
1138
 
      break;
1139
 
    case OP_STEP_MAX:
1140
 
      if (me == 0)
1141
 
    printf ("Testing GA_Step_max...");
1142
 
      if (type != C_DCPL || type != C_SCPL) {
1143
 
          /*NGA_Fill_patch (g_b, lo, hi, val2);*/
1144
 
          GA_Abs_value_patch (g_b, lo, hi);
1145
 
          GA_Step_max_patch (g_b, lo, hi, g_j, lo, hi, resultx);
1146
 
    /*
1147
 
    printf(" GA_Stepmax_patch type = %d, resultx = %le\n",type,resultx);
1148
 
    fflush(stdout);
1149
 
    */
1150
 
          /* 
1151
 
            It would be more robust to use GA_Elem_min_patch 
1152
 
            here to determine the minimum g_j value, but for
1153
 
            now we set it to -2.
1154
 
          */
1155
 
          aresulti = ((int)(GA_ABS(ival2)/GA_ABS(ival))) - resulti;
1156
 
          aresultd = GA_ABS(dval2/dval) - resultd;
1157
 
          aresultf = ((float)GA_ABS(fval2/fval)) - resultf;
1158
 
          aresultl = ((long)(GA_ABS(lval2)/GA_ABS(lval))) - resultl;
1159
 
      }
1160
 
      break;
1161
 
 
1162
 
    case OP_STEP_BOUND_INFO:
1163
 
      if (me == 0)
1164
 
    printf ("Testing GA_Step_bound_info...");
1165
 
      if (type != C_DCPL || type != C_SCPL) {
1166
 
          /*NGA_Fill_patch (g_b, lo, hi, val2);*/
1167
 
          GA_Abs_value_patch (g_b, lo, hi);
1168
 
          GA_Abs_value_patch (g_a, lo, hi);
1169
 
          /*GA_Abs_value_patch (g_j, lo, hi);*/
1170
 
          NGA_Fill_patch(g_c, lo, hi, val5);
1171
 
          GA_Step_bound_info_patch (g_b,lo,hi, g_j,lo,hi, g_a,lo,hi, g_c,lo,hi, boundminx,wolfeminx,boundmaxx);
1172
 
    /*
1173
 
    printf(" GA_Stepmax2_patch type = %d, resultx = %le\n",type,resultx);
1174
 
    fflush(stdout);
1175
 
    */
1176
 
          /* 
1177
 
          This is currently hardwired. would need to change if 
1178
 
          val, val2 or val5 change.
1179
 
          */
1180
 
    switch (type)
1181
 
      {
1182
 
      case C_INT:
1183
 
        awolfemini = ((int)(((int)1)/((int)2))) - wolfemini;
1184
 
        aboundmini = ((int)(((int)1)/((int)2))) - boundmini;
1185
 
        aboundmaxi = (int)GA_INFINITY_I - boundmaxi;
1186
 
        break;
1187
 
      case C_DBL:
1188
 
        awolfemind = (((double)1.0)/((double)2.0)) - wolfemind;
1189
 
        aboundmind = (((double)1.0)/((double)2.0)) - boundmind;
1190
 
        aboundmaxd = (double)GA_INFINITY_D - boundmaxd;
1191
 
        break;
1192
 
      case C_FLOAT:
1193
 
        awolfeminf = (((float)1.0)/((float)2.0))   - wolfeminf;
1194
 
        aboundminf = (((float)1.0)/((float)2.0))   - boundminf;
1195
 
        aboundmaxf = (float)GA_INFINITY_F - boundmaxf;
1196
 
        break;
1197
 
      case C_LONG:
1198
 
        awolfeminl = ((long)(((long)1)/((long)2))) - wolfeminl; 
1199
 
        aboundminl = ((long)(((long)1)/((long)2))) - boundminl; 
1200
 
        aboundmaxl = (long)GA_INFINITY_L - boundmaxl;
1201
 
        break;
1202
 
      default:
1203
 
        GA_Error ("GA_step_bound_info wrong data type.", type);
1204
 
      }
1205
 
      }
1206
 
      break;
1207
 
    default:
1208
 
      GA_Error ("test_function: wrong operation.", OP);
1209
 
    }
1210
 
  switch (type)
1211
 
    {
1212
 
    case C_INT:
1213
 
      alpha = &ai;
1214
 
      beta = &bi;
1215
 
      break;
1216
 
    case C_DCPL:
1217
 
      alpha = &adc;
1218
 
      beta = &bdc;
1219
 
      break;
1220
 
 
1221
 
    case C_SCPL:
1222
 
      alpha = &afc;
1223
 
      beta = &bfc;
1224
 
      break;
1225
 
 
1226
 
    case C_DBL:
1227
 
      alpha = &ad;
1228
 
      beta = &bd;
1229
 
      break;
1230
 
    case C_FLOAT:
1231
 
      alpha = &af;
1232
 
      beta = &bf;
1233
 
      break;
1234
 
    case C_LONG:
1235
 
      alpha = &al;
1236
 
      beta = &bl;
1237
 
      break;
1238
 
    default:
1239
 
      GA_Error ("wrong data type.", type);
1240
 
    }
1241
 
 
1242
 
  if (OP < 4) {
1243
 
    /* 
1244
 
      Binary operation. 
1245
 
    */
1246
 
    NGA_Add_patch (alpha, g_c, lo, hi, beta, g_d, lo, hi, g_e, lo, hi);
1247
 
    if (needs_scaled_result == 1) {
1248
 
      NGA_Add_patch (alpha, g_h, lo, hi, beta, g_i, lo, hi, g_j, lo, hi);
1249
 
      NGA_Add_patch (alpha, g_m, lo, hi, beta, g_n, lo, hi, g_n, lo, hi);
1250
 
    }
1251
 
  } else {
1252
 
 
1253
 
    /*
1254
 
      Unary operation.
1255
 
    */
1256
 
    if (OP < 7) {
1257
 
      NGA_Add_patch (alpha, g_a, lo, hi, beta, g_d, lo, hi, g_e, lo, hi);
1258
 
      if (needs_scaled_result == 1) {
1259
 
    NGA_Add_patch (alpha, g_f, lo, hi, beta, g_i, lo, hi, g_j, lo, hi);
1260
 
    NGA_Add_patch (alpha, g_k, lo, hi, beta, g_n, lo, hi, g_n, lo, hi);
1261
 
      }
1262
 
    }
1263
 
    /* 
1264
 
      Else it was a reduction operation (one of the step_max functions).
1265
 
    */
1266
 
  }
1267
 
 
1268
 
  switch (type)
1269
 
    {
1270
 
    case C_INT:
1271
 
      max = &imax;
1272
 
      break;
1273
 
    case C_DCPL:
1274
 
      max = &dcmax;
1275
 
      max2 = &dcmax2;
1276
 
      max3 = &dcmax3;
1277
 
      break;
1278
 
    case C_SCPL:
1279
 
      max = &fcmax;
1280
 
      max2 = &fcmax2;
1281
 
      max3 = &fcmax3;
1282
 
      break;
1283
 
    case C_DBL:
1284
 
      max = &dmax;
1285
 
      break;
1286
 
    case C_FLOAT:
1287
 
      max = &fmax;
1288
 
      break;
1289
 
    case C_LONG:
1290
 
      max = &lmax;
1291
 
      break;
1292
 
    default:
1293
 
      GA_Error ("wrong data type.", type);
1294
 
    }
1295
 
 
1296
 
  /*  
1297
 
    for unary and binary operators extract the maximum difference between
1298
 
    computed and correct solutions.
1299
 
  */
1300
 
  if (OP < 7) {
1301
 
    NGA_Select_elem (g_e, "max", max, index);
1302
 
    if (needs_scaled_result == 1) {
1303
 
      NGA_Select_elem (g_j, "max", max2, index2);
1304
 
      NGA_Select_elem (g_n, "max", max3, index3);
1305
 
    }
1306
 
  }
1307
 
  /*  NGA_Select_elem (g_e, "min", min, index);*/
1308
 
 
1309
 
  if (OP < 7) {
1310
 
    /* 
1311
 
      Binary or Unary operators.
1312
 
    */
1313
 
    switch (type)
1314
 
      {
1315
 
      double r, im;
1316
 
      float rf, imf;
1317
 
      case C_INT:
1318
 
      /*      result = (int)(imax - imin);*/
1319
 
      result = imax;
1320
 
      if (result != (int)0) result = 1;
1321
 
      break;
1322
 
      case C_DCPL:
1323
 
      /*
1324
 
      r = dcmax.real - dcmin.real;
1325
 
      im = dcmax.imag - dcmin.imag;
1326
 
      */
1327
 
      r = dcmax.real;
1328
 
      im = dcmax.imag;
1329
 
      if ((GA_ABS(r) + GA_ABS(im)) == (double)0.0) {
1330
 
       result = 0;
1331
 
      } else {
1332
 
       result = 1;
1333
 
      }
1334
 
      if (needs_scaled_result == 1) {
1335
 
      result2 = 0;
1336
 
       r = dcmax2.real;
1337
 
       im = dcmax2.imag;
1338
 
       if ((GA_ABS(r) + GA_ABS(im)) == (double)0.0) {
1339
 
         result2 = 0;
1340
 
       } else {
1341
 
         result2 = 1;
1342
 
       }
1343
 
       r = dcmax3.real;
1344
 
       im = dcmax3.imag;
1345
 
       if ((GA_ABS(r) + GA_ABS(im)) == (double)0.0) {
1346
 
         result3 = 0;
1347
 
       } else {
1348
 
         result3 = 1;
1349
 
       }
1350
 
      result = result | result2 | result3;
1351
 
      }
1352
 
      break;
1353
 
      case C_SCPL:
1354
 
      /*
1355
 
      rf = fcmax.real - fcmin.real;
1356
 
      imf = fcmax.imag - fcmin.imag;
1357
 
      */
1358
 
      rf = fcmax.real;
1359
 
      imf = fcmax.imag;
1360
 
      if ((GA_ABS(rf) + GA_ABS(imf)) == (float)0.0) {
1361
 
       result = 0;
1362
 
      } else {
1363
 
       result = 1;
1364
 
      }
1365
 
      if (needs_scaled_result == 1) {
1366
 
      result2 = 0;
1367
 
       rf = fcmax2.real;
1368
 
       imf = fcmax2.imag;
1369
 
       if ((GA_ABS(rf) + GA_ABS(imf)) == (float)0.0) {
1370
 
         result2 = 0;
1371
 
       } else {
1372
 
         result2 = 1;
1373
 
       }
1374
 
       rf = fcmax3.real;
1375
 
       imf = fcmax3.imag;
1376
 
       if ((GA_ABS(rf) + GA_ABS(imf)) == (float)0.0) {
1377
 
         result3 = 0;
1378
 
       } else {
1379
 
         result3 = 1;
1380
 
       }
1381
 
      result = result | result2 | result3;
1382
 
      }
1383
 
      break;
1384
 
      case C_DBL:
1385
 
      if (dmax == (double)0.0) {
1386
 
       result = 0;
1387
 
      } else {
1388
 
       result = 1;
1389
 
      }
1390
 
      break;
1391
 
      case C_FLOAT:
1392
 
      if (fmax == (float)0.0) {
1393
 
       result = 0;
1394
 
      } else {
1395
 
       result = 1;
1396
 
      }
1397
 
      break;
1398
 
      case C_LONG:
1399
 
      if (lmax == (long)0) {
1400
 
       result = 0;
1401
 
      } else {
1402
 
       result = 1;
1403
 
      }
1404
 
      break;
1405
 
      default:
1406
 
      GA_Error ("wrong data type.", type);
1407
 
      }
1408
 
  } else {
1409
 
    /*
1410
 
      A reduction operation, Step_max or Step_bound_info.
1411
 
    */
1412
 
    if (type == C_DCPL || type == C_SCPL) {
1413
 
      result = 0;
1414
 
    } else {
1415
 
      if (OP == OP_STEP_MAX) {
1416
 
    /* Step_max */
1417
 
    switch (type) 
1418
 
      {
1419
 
      case C_INT:
1420
 
        if (aresulti == 0) {
1421
 
          result = 0;
1422
 
        } else {
1423
 
          result = 1;
1424
 
        }
1425
 
        break;
1426
 
      case C_DBL:
1427
 
        if (aresultd == (double)0.0) {
1428
 
          result = 0;
1429
 
        } else {
1430
 
          result = 1;
1431
 
        }
1432
 
        break;
1433
 
      case C_FLOAT:
1434
 
        if (aresultf == (float)0.0) {
1435
 
          result = 0;
1436
 
        } else {
1437
 
          result = 1;
1438
 
        }
1439
 
        break;
1440
 
      case C_LONG:
1441
 
        if (aresultl == (long)0) {
1442
 
          result = 0;
1443
 
        } else {
1444
 
          result = 1;
1445
 
        }
1446
 
        break;
1447
 
      default:
1448
 
        GA_Error ("Stepmax op, wrong data type.", type);
1449
 
      }
1450
 
      } else {
1451
 
    /* OP = 8 so Step_bound_info */
1452
 
    switch (type) 
1453
 
      {
1454
 
      case C_INT:
1455
 
        if (awolfemini == 0) {
1456
 
          result = 0;
1457
 
        } else {
1458
 
          result = 1;
1459
 
        }
1460
 
        if (aboundmini == 0) {
1461
 
          result2 = 0;
1462
 
        } else {
1463
 
          result2 = 1;
1464
 
        }
1465
 
        if (aboundmaxi == 0) {
1466
 
          result3 = 0;
1467
 
        } else {
1468
 
          result3 = 1;
1469
 
        }
1470
 
        break;
1471
 
      case C_DBL:
1472
 
        if (awolfemind == ((double)0.0)) {
1473
 
          result = 0;
1474
 
        } else {
1475
 
          result = 1;
1476
 
        }
1477
 
        if (aboundmind == ((double)0.0)) {
1478
 
          result2 = 0;
1479
 
        } else {
1480
 
          result2 = 1;
1481
 
        }
1482
 
        if (aboundmaxd == ((double)0.0)) {
1483
 
          result3 = 0;
1484
 
        } else {
1485
 
          result3 = 1;
1486
 
        }
1487
 
        break;
1488
 
      case C_FLOAT:
1489
 
        if (awolfeminf == ((float)0.0)) {
1490
 
          result = 0;
1491
 
        } else {
1492
 
          result = 1;
1493
 
        }
1494
 
        if (aboundminf == ((float)0.0)) {
1495
 
          result2 = 0;
1496
 
        } else {
1497
 
          result2 = 1;
1498
 
        }
1499
 
        if (aboundmaxf == ((float)0.0)) {
1500
 
          result3 = 0;
1501
 
        } else {
1502
 
          result3 = 1;
1503
 
        }
1504
 
        break;
1505
 
      case C_LONG:
1506
 
        if (awolfeminl == ((long)0)) {
1507
 
          result = 0;
1508
 
        } else {
1509
 
          result = 1;
1510
 
        }
1511
 
        if (aboundminl == ((long)0)) {
1512
 
          result2 = 0;
1513
 
        } else {
1514
 
          result2 = 1;
1515
 
        }
1516
 
        if (aboundmaxl == ((long)0)) {
1517
 
          result3 = 0;
1518
 
        } else {
1519
 
          result3 = 1;
1520
 
        }
1521
 
        break;
1522
 
      default:
1523
 
        GA_Error ("Stepmax op, wrong data type.", type);
1524
 
      }
1525
 
    result = result | result2 | result3;
1526
 
      }
1527
 
    }
1528
 
  }
1529
 
  if (me == 0)
1530
 
    {
1531
 
      if (MISMATCHED (result, 0)) {
1532
 
    printf ("is not ok\n");
1533
 
    GA_Error("aborting", 1);
1534
 
      } else {
1535
 
    printf ("is ok.\n");
1536
 
      }
1537
 
    }
1538
 
 
1539
 
/*
1540
 
 NGA_Print_patch(g_a, lo, hi, 1);
1541
 
 NGA_Print_patch(g_d, lo, hi, 1);
1542
 
 NGA_Print_patch(g_e, lo, hi, 1);
1543
 
*/
1544
 
 
1545
 
  GA_Destroy (g_a);
1546
 
  GA_Destroy (g_b);
1547
 
  GA_Destroy (g_c);
1548
 
  GA_Destroy (g_d);
1549
 
  GA_Destroy (g_e);
1550
 
  GA_Destroy (g_f);
1551
 
  GA_Destroy (g_g);
1552
 
  GA_Destroy (g_h);
1553
 
  GA_Destroy (g_i);
1554
 
  GA_Destroy (g_j);
1555
 
  GA_Destroy (g_k);
1556
 
  GA_Destroy (g_l);
1557
 
  GA_Destroy (g_m);
1558
 
  GA_Destroy (g_n);
1559
 
 
1560
 
  return ok;
1561
 
}
1562
 
 
1563
 
int
1564
 
main (argc, argv)
1565
 
     int argc;
1566
 
     char **argv;
1567
 
{
1568
 
  int heap = 20000, stack = 20000;
1569
 
  int me, nproc;
1570
 
  int d, op;
1571
 
  int ok = 1;
1572
 
 
1573
 
  MP_INIT(argc,argv);
1574
 
  GA_INIT(argc,argv);        /* initialize GA */
1575
 
  me = GA_Nodeid ();
1576
 
  nproc = GA_Nnodes ();
1577
 
  if (me == 0)
1578
 
    {
1579
 
      if (GA_Uses_fapi ())
1580
 
    GA_Error ("Program runs with C array API only", 1);
1581
 
      printf ("Using %ld processes\n", (long) nproc);
1582
 
      fflush (stdout);
1583
 
    }
1584
 
 
1585
 
  heap /= nproc;
1586
 
  stack /= nproc;
1587
 
  if (!MA_init (C_DBL, stack, heap))
1588
 
    GA_Error ("MA_init failed", stack + heap);    /* initialize memory allocator */
1589
 
 
1590
 
 
1591
 
  /* op = 8;*/
1592
 
  for (op = 0; op < 9; op++)
1593
 
    {
1594
 
      /*for (d = 1; d < 2; d++)*/
1595
 
      for (d = 1; d < 4; d++)
1596
 
    {
1597
 
      if (me == 0) 
1598
 
        printf ("\n\ndim =%d\n\n", d);
1599
 
      if (me == 0) 
1600
 
        printf ("\ndata type: INT\t\t");
1601
 
      ok = test_fun (C_INT, d, op);
1602
 
      if (me == 0) 
1603
 
        printf ("\ndata type: double\t");
1604
 
      ok = test_fun (C_DBL, d, op);
1605
 
      if (me == 0)
1606
 
        printf ("\ndata type: float\t");
1607
 
      ok = test_fun (C_FLOAT, d, op);
1608
 
      if (me == 0)
1609
 
        printf ("\ndata type: long\t\t");
1610
 
      ok = test_fun (C_LONG, d, op);
1611
 
      if (op < 7) {
1612
 
        if (me == 0)
1613
 
          printf ("\ndata type: double complex\t");
1614
 
        ok = test_fun (C_DCPL, d, op);
1615
 
      }
1616
 
    }
1617
 
    }
1618
 
 
1619
 
  printf("\nAll tests successful\n");
1620
 
 
1621
 
  GA_Terminate();
1622
 
  
1623
 
  MP_FINALIZE();
1624
 
 
1625
 
  return 0;
1626
 
}
1627
 
 
1628
 
/*\ FILL IN ARRAY WITH Varying VALUEs. (from 0 to product of dims-1).
1629
 
    For complex arrays make the real and imaginary parts equal.
1630
 
\*/
1631
 
void nga_vfill_patch(Integer *g_a, Integer *lo, Integer *hi)
1632
 
{
1633
 
    Integer i, j;
1634
 
    Integer ndim, dims[MAXDIM], type;
1635
 
    Integer loA[MAXDIM], hiA[MAXDIM], ld[MAXDIM];
1636
 
    void *data_ptr;
1637
 
    Integer idx, n1dim;
1638
 
    Integer bvalue[MAXDIM], bunit[MAXDIM], baseld[MAXDIM];
1639
 
    Integer me= pnga_nodeid();
1640
 
    int local_sync_begin,local_sync_end;
1641
 
   
1642
 
    local_sync_begin = _ga_sync_begin; local_sync_end = _ga_sync_end;
1643
 
    _ga_sync_begin = 1; _ga_sync_end=1; /*remove any previous masking*/
1644
 
    if(local_sync_begin)GA_Sync(); 
1645
 
 
1646
 
    GA_PUSH_NAME("nga_vfill_patch");
1647
 
    
1648
 
    pnga_inquire(*g_a,  &type, &ndim, dims);
1649
 
 
1650
 
    /* get limits of VISIBLE patch */ 
1651
 
    pnga_distribution(*g_a, me, loA, hiA);
1652
 
    
1653
 
    /*  determine subset of my local patch to access  */
1654
 
    /*  Output is in loA and hiA */
1655
 
    if(pnga_patch_intersect(lo, hi, loA, hiA, ndim)){
1656
 
 
1657
 
        /* get data_ptr to corner of patch */
1658
 
        /* ld are leading dimensions INCLUDING ghost cells */
1659
 
        pnga_access_ptr(*g_a, loA, hiA, &data_ptr, ld);
1660
 
 
1661
 
        /* number of n-element of the first dimension */
1662
 
        n1dim = 1; for(i=1; i<ndim; i++) n1dim *= (hiA[i] - loA[i] + 1);
1663
 
        
1664
 
        /* calculate the destination indices */
1665
 
        bvalue[0] = 0; bvalue[1] = 0; bunit[0] = 1; bunit[1] = 1;
1666
 
        /* baseld[0] = ld[0]
1667
 
         * baseld[1] = ld[0] * ld[1]
1668
 
         * baseld[2] = ld[0] * ld[1] * ld[2] .....
1669
 
         */
1670
 
        baseld[0] = ld[0]; baseld[1] = baseld[0] *ld[1];
1671
 
        for(i=2; i<ndim; i++) {
1672
 
            bvalue[i] = 0;
1673
 
            bunit[i] = bunit[i-1] * (hiA[i-1] - loA[i-1] + 1);
1674
 
            baseld[i] = baseld[i-1] * ld[i];
1675
 
        }
1676
 
 
1677
 
        switch (type){
1678
 
            case C_INT:
1679
 
                for(i=0; i<n1dim; i++) {
1680
 
                    idx = 0;
1681
 
                    for(j=1; j<ndim; j++) {
1682
 
                        idx += bvalue[j] * baseld[j-1];
1683
 
                        if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1684
 
                        if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1685
 
                    }
1686
 
                    for(j=0; j<(hiA[0]-loA[0]+1); j++)
1687
 
                        ((int *)data_ptr)[idx+j] = (int)(idx+j);
1688
 
                }
1689
 
                break;
1690
 
            case C_DCPL:
1691
 
                for(i=0; i<n1dim; i++) {
1692
 
                    idx = 0;
1693
 
                    for(j=1; j<ndim; j++) {
1694
 
                        idx += bvalue[j] * baseld[j-1];
1695
 
                        if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1696
 
                        if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1697
 
                    }
1698
 
                    
1699
 
                    for(j=0; j<(hiA[0]-loA[0]+1); j++) {
1700
 
                        ((DoubleComplex *)data_ptr)[idx+j].real = (double)(idx+j);
1701
 
                        ((DoubleComplex *)data_ptr)[idx+j].imag = (double)(idx+j);
1702
 
                    }
1703
 
                }
1704
 
                
1705
 
                break;
1706
 
            case C_SCPL:
1707
 
                for(i=0; i<n1dim; i++) {
1708
 
                    idx = 0;
1709
 
                    for(j=1; j<ndim; j++) {
1710
 
                        idx += bvalue[j] * baseld[j-1];
1711
 
                        if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1712
 
                        if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1713
 
                    }
1714
 
                    
1715
 
                    for(j=0; j<(hiA[0]-loA[0]+1); j++) {
1716
 
                        ((SingleComplex *)data_ptr)[idx+j].real = (float)(idx+j);
1717
 
                        ((SingleComplex *)data_ptr)[idx+j].imag = (float)(idx+j);
1718
 
                    }
1719
 
                }
1720
 
                
1721
 
                break;
1722
 
            case C_DBL:
1723
 
                for(i=0; i<n1dim; i++) {
1724
 
                    idx = 0;
1725
 
                    for(j=1; j<ndim; j++) {
1726
 
                        idx += bvalue[j] * baseld[j-1];
1727
 
                        if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1728
 
                        if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1729
 
                    }
1730
 
                    
1731
 
                    for(j=0; j<(hiA[0]-loA[0]+1); j++) 
1732
 
              ((double*)data_ptr)[idx+j] = (double)(idx+j);
1733
 
                }
1734
 
                break;
1735
 
            case C_FLOAT:
1736
 
                for(i=0; i<n1dim; i++) {
1737
 
                    idx = 0;
1738
 
                    for(j=1; j<ndim; j++) {
1739
 
                        idx += bvalue[j] * baseld[j-1];
1740
 
                        if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1741
 
                        if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1742
 
                    }
1743
 
                    for(j=0; j<(hiA[0]-loA[0]+1); j++)
1744
 
                        ((float *)data_ptr)[idx+j] = (float)(idx+j);
1745
 
                }
1746
 
                break;     
1747
 
            case C_LONG:
1748
 
                for(i=0; i<n1dim; i++) {
1749
 
                    idx = 0;
1750
 
                    for(j=1; j<ndim; j++) {
1751
 
                        idx += bvalue[j] * baseld[j-1];
1752
 
                        if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1753
 
                        if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1754
 
                    }
1755
 
                    for(j=0; j<(hiA[0]-loA[0]+1); j++)
1756
 
                        ((long *)data_ptr)[idx+j] = (long)(idx+j);
1757
 
                } 
1758
 
                break;                          
1759
 
            default: GA_Error(" wrong data type ",type);
1760
 
        }
1761
 
        
1762
 
        /* release access to the data */
1763
 
        pnga_release_update(*g_a, loA, hiA);
1764
 
    }
1765
 
    GA_POP_NAME;
1766
 
    if(local_sync_end)GA_Sync();
1767
 
}
1768
 
/*\ Utility function to actually set positive/negative values
1769
 
\*/
1770
 
void ngai_do_pnfill_patch(Integer type, Integer ndim, Integer *loA, Integer *hiA,
1771
 
                          Integer *ld, void *data_ptr)
1772
 
{
1773
 
  Integer i, j;
1774
 
  Integer idx, n1dim;
1775
 
  Integer bvalue[MAXDIM], bunit[MAXDIM], baseld[MAXDIM];
1776
 
  /* number of n-element of the first dimension */
1777
 
  n1dim = 1; for(i=1; i<ndim; i++) n1dim *= (hiA[i] - loA[i] + 1);
1778
 
 
1779
 
  /* calculate the destination indices */
1780
 
  bvalue[0] = 0; bvalue[1] = 0; bunit[0] = 1; bunit[1] = 1;
1781
 
  /* baseld[0] = ld[0]
1782
 
   * baseld[1] = ld[0] * ld[1]
1783
 
   * baseld[2] = ld[0] * ld[1] * ld[2] .....
1784
 
   */
1785
 
  baseld[0] = ld[0]; baseld[1] = baseld[0] *ld[1];
1786
 
  for(i=2; i<ndim; i++) {
1787
 
    bvalue[i] = 0;
1788
 
    bunit[i] = bunit[i-1] * (hiA[i-1] - loA[i-1] + 1);
1789
 
    baseld[i] = baseld[i-1] * ld[i];
1790
 
  }
1791
 
 
1792
 
  switch (type){
1793
 
    case C_INT:
1794
 
      for(i=0; i<n1dim; i++) {
1795
 
        idx = 0;
1796
 
        for(j=1; j<ndim; j++) {
1797
 
          idx += bvalue[j] * baseld[j-1];
1798
 
          if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1799
 
          if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1800
 
        }
1801
 
 
1802
 
        for(j=0; j<(hiA[0]-loA[0]+1); j++)
1803
 
          ((int *)data_ptr)[idx+j] = (int)(((idx+j)&3)-2);
1804
 
      }
1805
 
      break;
1806
 
    case C_DCPL:
1807
 
      for(i=0; i<n1dim; i++) {
1808
 
        idx = 0;
1809
 
        for(j=1; j<ndim; j++) {
1810
 
          idx += bvalue[j] * baseld[j-1];
1811
 
          if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1812
 
          if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1813
 
        }
1814
 
 
1815
 
        for(j=0; j<(hiA[0]-loA[0]+1); j++) {
1816
 
          ((DoubleComplex *)data_ptr)[idx+j].real = (double)(((idx+j)&3)-2);
1817
 
          ((DoubleComplex *)data_ptr)[idx+j].imag = (double)(((idx+j)&3)-2);
1818
 
        }
1819
 
      }
1820
 
      break;
1821
 
    case C_SCPL:
1822
 
      for(i=0; i<n1dim; i++) {
1823
 
        idx = 0;
1824
 
        for(j=1; j<ndim; j++) {
1825
 
          idx += bvalue[j] * baseld[j-1];
1826
 
          if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1827
 
          if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1828
 
        }
1829
 
 
1830
 
        for(j=0; j<(hiA[0]-loA[0]+1); j++) {
1831
 
          ((SingleComplex *)data_ptr)[idx+j].real = (float)(((idx+j)&3)-2);
1832
 
          ((SingleComplex *)data_ptr)[idx+j].imag = (float)(((idx+j)&3)-2);
1833
 
        }
1834
 
      }
1835
 
      break;
1836
 
    case C_DBL:
1837
 
      for(i=0; i<n1dim; i++) {
1838
 
        idx = 0;
1839
 
        for(j=1; j<ndim; j++) {
1840
 
          idx += bvalue[j] * baseld[j-1];
1841
 
          if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1842
 
          if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1843
 
        }
1844
 
 
1845
 
        for(j=0; j<(hiA[0]-loA[0]+1); j++) 
1846
 
          ((double*)data_ptr)[idx+j] = (double)(((idx+j)&3)-2);
1847
 
      }
1848
 
      break;
1849
 
    case C_FLOAT:
1850
 
      for(i=0; i<n1dim; i++) {
1851
 
        idx = 0;
1852
 
        for(j=1; j<ndim; j++) {
1853
 
          idx += bvalue[j] * baseld[j-1];
1854
 
          if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1855
 
          if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1856
 
        }
1857
 
        for(j=0; j<(hiA[0]-loA[0]+1); j++)
1858
 
          ((float *)data_ptr)[idx+j] = (float)(((idx+j)&3)-2);
1859
 
      }
1860
 
      break;     
1861
 
    case C_LONG:
1862
 
      for(i=0; i<n1dim; i++) {
1863
 
        idx = 0;
1864
 
        for(j=1; j<ndim; j++) {
1865
 
          idx += bvalue[j] * baseld[j-1];
1866
 
          if(((i+1) % bunit[j]) == 0) bvalue[j]++;
1867
 
          if(bvalue[j] > (hiA[j]-loA[j])) bvalue[j] = 0;
1868
 
        }
1869
 
        for(j=0; j<(hiA[0]-loA[0]+1); j++)
1870
 
          ((long *)data_ptr)[idx+j] = (long)(((idx+j)&3)-2);
1871
 
      } 
1872
 
      break;                          
1873
 
    default: GA_Error(" wrong data type ",type);
1874
 
  }
1875
 
 
1876
 
}
1877
 
 
1878
 
/*\ FILL IN ARRAY WITH Varying positive and negative VALUEs. 
1879
 
    (from -2 to 1).
1880
 
    For complex arrays make the real and imaginary parts equal.
1881
 
\*/
1882
 
void nga_pnfill_patch(Integer *g_a, Integer *lo, Integer *hi)
1883
 
{
1884
 
  Integer i;
1885
 
  Integer ndim, dims[MAXDIM], type;
1886
 
  Integer loA[MAXDIM], hiA[MAXDIM], ld[MAXDIM];
1887
 
  void *data_ptr;
1888
 
  Integer me= pnga_nodeid();
1889
 
  Integer num_blocks;
1890
 
  int local_sync_begin,local_sync_end;
1891
 
 
1892
 
  local_sync_begin = _ga_sync_begin; local_sync_end = _ga_sync_end;
1893
 
  _ga_sync_begin = 1; _ga_sync_end=1; /*remove any previous masking*/
1894
 
  if(local_sync_begin)GA_Sync(); 
1895
 
 
1896
 
  GA_PUSH_NAME("nga_pnfill_patch");
1897
 
 
1898
 
  pnga_inquire(*g_a,  &type, &ndim, dims);
1899
 
 
1900
 
  num_blocks = pnga_total_blocks(*g_a);
1901
 
 
1902
 
  if (num_blocks < 0) {
1903
 
    /* get limits of VISIBLE patch */ 
1904
 
    pnga_distribution(*g_a, me, loA, hiA);
1905
 
 
1906
 
    /*  determine subset of my local patch to access  */
1907
 
    /*  Output is in loA and hiA */
1908
 
    if(pnga_patch_intersect(lo, hi, loA, hiA, ndim)){
1909
 
 
1910
 
      /* get data_ptr to corner of patch */
1911
 
      /* ld are leading dimensions INCLUDING ghost cells */
1912
 
      pnga_access_ptr(*g_a, loA, hiA, &data_ptr, ld);
1913
 
 
1914
 
      ngai_do_pnfill_patch(type, ndim, loA, hiA, ld, data_ptr);
1915
 
 
1916
 
      /* release access to the data */
1917
 
      pnga_release_update(*g_a, loA, hiA);
1918
 
    }
1919
 
  } else {
1920
 
    Integer offset, j, jtmp, chk;
1921
 
    Integer loS[MAXDIM];
1922
 
    Integer nproc = pnga_nnodes();
1923
 
    /* using simple block-cyclic data distribution */
1924
 
    if (!pnga_uses_proc_grid(*g_a)){
1925
 
      for (i=me; i<num_blocks; i += nproc) {
1926
 
        /* get limits of patch */
1927
 
        pnga_distribution(*g_a, i, loA, hiA);
1928
 
 
1929
 
        /* loA is changed by pnga_patch_intersect, so
1930
 
           save a copy */
1931
 
        for (j=0; j<ndim; j++) {
1932
 
          loS[j] = loA[j];
1933
 
        }
1934
 
 
1935
 
        /*  determine subset of my local patch to access  */
1936
 
        /*  Output is in loA and hiA */
1937
 
        if(pnga_patch_intersect(lo, hi, loA, hiA, ndim)){
1938
 
 
1939
 
          /* get data_ptr to corner of patch */
1940
 
          /* ld are leading dimensions for block */
1941
 
          pnga_access_block_ptr(*g_a, i, &data_ptr, ld);
1942
 
 
1943
 
          /* Check for partial overlap */
1944
 
          chk = 1;
1945
 
          for (j=0; j<ndim; j++) {
1946
 
            if (loS[j] < loA[j]) {
1947
 
              chk=0;
1948
 
              break;
1949
 
            }
1950
 
          }
1951
 
          if (!chk) {
1952
 
            /* Evaluate additional offset for pointer */
1953
 
            offset = 0;
1954
 
            jtmp = 1;
1955
 
            for (j=0; j<ndim-1; j++) {
1956
 
              offset += (loA[j]-loS[j])*jtmp;
1957
 
              jtmp *= ld[j];
1958
 
            }
1959
 
            offset += (loA[ndim-1]-loS[ndim-1])*jtmp;
1960
 
            switch (type){
1961
 
              case C_INT:
1962
 
                data_ptr = (void*)((int*)data_ptr + offset);
1963
 
                break;
1964
 
              case C_DCPL:
1965
 
                data_ptr = (void*)((double*)data_ptr + 2*offset);
1966
 
                break;
1967
 
              case C_SCPL:
1968
 
                data_ptr = (void*)((float*)data_ptr + 2*offset);
1969
 
                break;
1970
 
              case C_DBL:
1971
 
                data_ptr = (void*)((double*)data_ptr + offset);
1972
 
                break;
1973
 
              case C_FLOAT:
1974
 
                data_ptr = (void*)((float*)data_ptr + offset);
1975
 
                break;
1976
 
              case C_LONG:
1977
 
                data_ptr = (void*)((long*)data_ptr + offset);
1978
 
                break;
1979
 
              default: GA_Error(" wrong data type ",type);
1980
 
            }
1981
 
          }
1982
 
          /* fill in patch */
1983
 
          ngai_do_pnfill_patch(type, ndim, loA, hiA, ld, data_ptr);
1984
 
 
1985
 
          /* release access to the data */
1986
 
          pnga_release_update_block(*g_a, i);
1987
 
        }
1988
 
      }
1989
 
    } else {
1990
 
      /* using scalapack block-cyclic data distribution */
1991
 
      Integer proc_index[MAXDIM], index[MAXDIM];
1992
 
      Integer topology[MAXDIM];
1993
 
      Integer blocks[MAXDIM], block_dims[MAXDIM];
1994
 
      pnga_get_proc_index(*g_a, me, proc_index);
1995
 
      pnga_get_proc_index(*g_a, me, index);
1996
 
      pnga_get_block_info(*g_a, blocks, block_dims);
1997
 
      pnga_get_proc_grid(*g_a, topology);
1998
 
      while (index[ndim-1] < blocks[ndim-1]) {
1999
 
        /* find bounding coordinates of block */
2000
 
        chk = 1;
2001
 
        for (i = 0; i < ndim; i++) {
2002
 
          loA[i] = index[i]*block_dims[i]+1;
2003
 
          hiA[i] = (index[i] + 1)*block_dims[i];
2004
 
          if (hiA[i] > dims[i]) hiA[i] = dims[i];
2005
 
          if (hiA[i] < loA[i]) chk = 0;
2006
 
        }
2007
 
 
2008
 
        /* loA is changed by pnga_patch_intersect, so
2009
 
           save a copy */
2010
 
        for (j=0; j<ndim; j++) {
2011
 
          loS[j] = loA[j];
2012
 
        }
2013
 
 
2014
 
        /*  determine subset of my local patch to access  */
2015
 
        /*  Output is in loA and hiA */
2016
 
        if(pnga_patch_intersect(lo, hi, loA, hiA, ndim)){
2017
 
 
2018
 
          /* get data_ptr to corner of patch */
2019
 
          /* ld are leading dimensions for block */
2020
 
          pnga_access_block_grid_ptr(*g_a, index, &data_ptr, ld);
2021
 
 
2022
 
          /* Check for partial overlap */
2023
 
          chk = 1;
2024
 
          for (j=0; j<ndim; j++) {
2025
 
            if (loS[j] < loA[j]) {
2026
 
              chk=0;
2027
 
              break;
2028
 
            }
2029
 
          }
2030
 
          if (!chk) {
2031
 
            /* Evaluate additional offset for pointer */
2032
 
            offset = 0;
2033
 
            jtmp = 1;
2034
 
            for (j=0; j<ndim-1; j++) {
2035
 
              offset += (loA[j]-loS[j])*jtmp;
2036
 
              jtmp *= ld[j];
2037
 
            }
2038
 
            offset += (loA[ndim-1]-loS[ndim-1])*jtmp;
2039
 
            switch (type){
2040
 
              case C_INT:
2041
 
                data_ptr = (void*)((int*)data_ptr + offset);
2042
 
                break;
2043
 
              case C_DCPL:
2044
 
                data_ptr = (void*)((double*)data_ptr + 2*offset);
2045
 
                break;
2046
 
              case C_SCPL:
2047
 
                data_ptr = (void*)((float*)data_ptr + 2*offset);
2048
 
                break;
2049
 
              case C_DBL:
2050
 
                data_ptr = (void*)((double*)data_ptr + offset);
2051
 
                break;
2052
 
              case C_FLOAT:
2053
 
                data_ptr = (void*)((float*)data_ptr + offset);
2054
 
                break;
2055
 
              case C_LONG:
2056
 
                data_ptr = (void*)((long*)data_ptr + offset);
2057
 
                break;
2058
 
              default: GA_Error(" wrong data type ",type);
2059
 
            }
2060
 
          }
2061
 
          /* fill in patch */
2062
 
          ngai_do_pnfill_patch(type, ndim, loA, hiA, ld, data_ptr);
2063
 
 
2064
 
          /* release access to the data */
2065
 
          pnga_release_update_block_grid(*g_a, index);
2066
 
        }
2067
 
        /* increment index to get next block on processor */
2068
 
        index[0] += topology[0];
2069
 
        for (i = 0; i < ndim; i++) {
2070
 
          if (index[i] >= blocks[i] && i<ndim-1) {
2071
 
            index[i] = proc_index[i];
2072
 
            index[i+1] += topology[i+1];
2073
 
          }
2074
 
        }
2075
 
      }
2076
 
    }
2077
 
  }
2078
 
  GA_POP_NAME;
2079
 
  if(local_sync_end)GA_Sync();
2080
 
}
2081
 
 
2082