← 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-2/global/testing/ga_lu.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-2/global/testing/ga_lu.c
 
1
#if HAVE_CONFIG_H
 
2
#   include "config.h"
 
3
#endif
 
4
 
 
5
 
 
6
#if HAVE_STDIO_H
 
7
#   include <stdio.h>
 
8
#endif
 
9
#if HAVE_STDLIB_H
 
10
#   include <stdlib.h>
 
11
#endif
 
12
 
 
13
#include "ga.h"
 
14
#include "macdecls.h"
 
15
#include "mp3.h"
 
16
#include "galinalg.h"
 
17
 
 
18
#define BLOCK_CYCLIC
 
19
#define BLOCK_SIZE 500
 
20
/*#define USE_SCALAPACK_DISTR*/
 
21
 
 
22
#define DEBUG 0
 
23
 
 
24
static int nprocs, me;
 
25
 
 
26
static void init_array(double *a, int n) 
 
27
{
 
28
    int i, j;
 
29
    double max=0.0; /* for pivoting */
 
30
    
 
31
    for(i=0; i<n; i++) 
 
32
    {
 
33
       for(j=0; j<n; j++) 
 
34
       {
 
35
          a[i*n+j] = (double)rand()/RAND_MAX; /* (i*n+j+1); */
 
36
          if(max<a[i*n+j]) max = a[i*n+j];
 
37
       }
 
38
    }
 
39
 
 
40
    /* Forcing the diagonal is always max to avoid pivoting*/
 
41
    for(i=0; i<n; i++) 
 
42
    {
 
43
       a[i*n+i] = max + (double)i;
 
44
    }
 
45
}
 
46
 
 
47
void copy_array(double *a, double *b, int n) 
 
48
{
 
49
    int i, j;
 
50
    
 
51
    for(i=0; i<n; i++) 
 
52
    {
 
53
       for(j=0; j<n; j++) 
 
54
       {
 
55
          b[i*n+j] = a[i*n+j]; 
 
56
       }
 
57
    }
 
58
}
 
59
 
 
60
void print_array(double *a, int n) 
 
61
{
 
62
    int i, j;
 
63
    
 
64
    printf("Print Matrix of size %d x %d:\n", n, n);
 
65
    for(i=0; i<n; i++) 
 
66
    {
 
67
       for(j=0; j<n; j++) 
 
68
       {
 
69
          if(a[i*n+j]<0) printf("%2.4e ",  a[i*n+j]);
 
70
          else           printf(" %2.4e ", a[i*n+j]);
 
71
       }
 
72
       printf("\n");
 
73
    }
 
74
    printf("\n");
 
75
}
 
76
 
 
77
/* print U, i.e. upper triangular */
 
78
void print_U(double *a, int n) 
 
79
{
 
80
    int i, j;
 
81
    
 
82
    printf("Print Upper triangular matrix of size %d x %d:\n", n, n);
 
83
    for(i=0; i<n; i++) 
 
84
    {
 
85
       for(j=0; j<i; j++) 
 
86
       {
 
87
          printf(" %2.4e ", 0.0);
 
88
       }
 
89
       
 
90
       for(j=i; j<n; j++) 
 
91
       {
 
92
          if(a[i*n+j]<0) printf("%2.4e ",  a[i*n+j]);
 
93
          else           printf(" %2.4e ", a[i*n+j]);
 
94
       }
 
95
       
 
96
       printf("\n");
 
97
    }
 
98
    printf("\n");
 
99
}
 
100
 
 
101
/* print L, i.e. unit lower triangular */
 
102
void print_L(double *a, int n) 
 
103
{
 
104
    int i, j;
 
105
    
 
106
    printf("Print Unit Lower triangular matrix of size %d x %d:\n", n, n);
 
107
    for(i=0; i<n; i++) 
 
108
    {
 
109
       for(j=0; j<i; j++) 
 
110
       {
 
111
          if(a[i*n+j]<0) printf("%2.4e ",  a[i*n+j]);
 
112
          else           printf(" %2.4e ", a[i*n+j]);
 
113
       }
 
114
       for(j=i; j<n; j++) 
 
115
       {
 
116
          if(i==j) printf(" %2.4e ", 1.0);
 
117
          else     printf(" %2.4e ", 0.0);
 
118
       }
 
119
 
 
120
       printf("\n");
 
121
    }
 
122
    printf("\n");
 
123
}
 
124
 
 
125
void lu_basic(double *a, int n) {
 
126
 
 
127
    int i, j, k;
 
128
    
 
129
    for (i=0; i<n-1; i++) 
 
130
    {
 
131
       /* scale lower matrix column */
 
132
       for (j=i+1; j<n; j++) 
 
133
       {
 
134
          a[j*n+i] = a[j*n+i]/a[i*n+i]; /* a[j][i] = a[j][i]/a[i][i] */
 
135
       }
 
136
       
 
137
       /* rank-one update of matrix */
 
138
       for (j=i+1; j<n; j++) 
 
139
       {
 
140
          for (k=i+1; k<n; k++) 
 
141
          {
 
142
             /*a[k][j] = a[k][j] - a[i][j]/a[k][i] */
 
143
             a[k*n+j]  = a[k*n+j] - a[i*n+j]/a[k*n+i];
 
144
          }   
 
145
       }   
 
146
    }
 
147
}
 
148
 
 
149
int lu_lapack(double *a, int n) 
 
150
{
 
151
    int i, j;
 
152
    double *aa=NULL;
 
153
    BlasInt *ipiv=NULL, info;
 
154
    BlasInt ld=(BlasInt)n;
 
155
    BlasInt N=(BlasInt)n;
 
156
    
 
157
    aa = (double*)malloc(n*n*sizeof(double));
 
158
    ipiv = (BlasInt*)malloc(n*sizeof(BlasInt));
 
159
 
 
160
    /* row-major to column-major (NOTE: dgetrf_ is a fortran function) */
 
161
    for(i=0; i<n; i++) 
 
162
    {
 
163
       for(j=0; j<n; j++) 
 
164
       {
 
165
          aa[j*n+i] = a[i*n+j];
 
166
       }
 
167
    }
 
168
    
 
169
    LAPACK_DGETRF(&N, &N, aa, &ld, ipiv, &info); /* LAPACK's LU */
 
170
    
 
171
    /* column-major to row-major */
 
172
    for(i=0; i<n; i++) 
 
173
    {
 
174
       for(j=0; j<n; j++) 
 
175
       {
 
176
          a[i*n+j] = aa[j*n+i];
 
177
       }
 
178
    }
 
179
    
 
180
#if DEBUG
 
181
    printf("ipiv[] = [ ");
 
182
    for(i=0; i<n; i++) printf("%ld ", ipiv[i]);
 
183
    printf("]\n");
 
184
#endif
 
185
    
 
186
    free(aa);
 
187
    return (int)info;
 
188
}
 
189
 
 
190
void lu(double *A, int matrix_size) 
 
191
{
 
192
    double *a=NULL, *a_verify=NULL;
 
193
    int info;
 
194
    
 
195
    a        = (double*)malloc(matrix_size*matrix_size*sizeof(double));
 
196
    a_verify = (double*)malloc(matrix_size*matrix_size*sizeof(double));
 
197
    if(a == NULL || a_verify == NULL) 
 
198
    {
 
199
       GA_Error("lu(): malloc failed", matrix_size*matrix_size*sizeof(double));
 
200
    }
 
201
    
 
202
    copy_array(A, a, matrix_size);
 
203
    copy_array(A, a_verify, matrix_size);
 
204
    
 
205
#if 0
 
206
    printf("\nDoing LU Factorization\n\n");
 
207
    lu_basic(a, matrix_size);
 
208
    printf("LU = \n");
 
209
    print_array(a, matrix_size);
 
210
    
 
211
    printf("\nDoing LAPACK's LU Factorization\n\n");
 
212
#endif
 
213
    info = lu_lapack(a_verify, matrix_size);
 
214
 
 
215
#if DEBUG
 
216
    printf("LU = ");
 
217
    print_array(a_verify, matrix_size);
 
218
#endif
 
219
    
 
220
    if(info!=0) 
 
221
    {
 
222
       printf("\nError: dgetrf() of lapack is NOT successful (INFO=%d)\n\n",
 
223
              info);
 
224
       printf("NOTE:\n INFO=0:  successful exit\n INFO<0:  if INFO = -i, the i-th argument had an illegal value\n INFO>0:  if INFO = i, U(i,i) is exactly zero. The factorization\n has been completed, but the factor U is exactly singular,  and\n division by zero will occur if it is used to solve a system of\n equations.\n");
 
225
       exit(0);
 
226
    }
 
227
    
 
228
 
 
229
    free(a_verify);
 
230
    free(a);
 
231
 
 
232
    
 
233
}
 
234
 
 
235
void dtrsm_lapack(double *a, double *b, int n,
 
236
                  char side, char uplo, char transa, char diag) 
 
237
{
 
238
    int i, j;
 
239
    double *aa=NULL;
 
240
    double *bb=NULL;
 
241
    BlasInt ld = (BlasInt)n;
 
242
    BlasInt N  = (BlasInt)n;
 
243
 
 
244
    aa = (double*)malloc(n*n*sizeof(double));
 
245
    bb = (double*)malloc(n*n*sizeof(double));
 
246
    
 
247
    /* row-major to column-major (NOTE: dgetrf_ is a fortran function) */
 
248
    for(i=0; i<n; i++) 
 
249
    {
 
250
       for(j=0; j<n; j++) 
 
251
       {
 
252
          aa[j*n+i] = a[i*n+j];
 
253
          bb[j*n+i] = b[i*n+j];
 
254
       }
 
255
    }
 
256
    
 
257
    {
 
258
       double alpha= 1.0;
 
259
       
 
260
       LAPACK_DTRSM(&side, &uplo, &transa, &diag, &N, &N, &alpha, aa, &ld, bb, &ld);
 
261
    }   
 
262
    
 
263
    /* column-major to row-major */
 
264
    for(i=0; i<n; i++) 
 
265
    {
 
266
       for(j=0; j<n; j++) 
 
267
       {
 
268
          a[i*n+j] = aa[j*n+i];
 
269
          b[i*n+j] = bb[j*n+i];
 
270
       }
 
271
    }
 
272
    
 
273
    free(aa);
 
274
    free(bb);    
 
275
}
 
276
 
 
277
void transpose(double *a, int n) 
 
278
{
 
279
    int i, j;
 
280
    double *aa=NULL;
 
281
    
 
282
    aa = (double*)malloc(n*n*sizeof(double));
 
283
    copy_array(a, aa, n);
 
284
    
 
285
    /* row-major to column-major (NOTE: dgetrf_ is a fortran function) */
 
286
    for(i=0; i<n; i++) 
 
287
    {
 
288
       for(j=0; j<n; j++) 
 
289
       {
 
290
          a[j*n+i] = aa[i*n+j];
 
291
       }
 
292
    }
 
293
}
 
294
 
 
295
 
 
296
/* input is matrix size */
 
297
void ga_lu(double *A, int matrix_size) 
 
298
{
 
299
    int g_a, g_b, dims[2], type=C_DBL;
 
300
    int lo[2], hi[2], ld;
 
301
    int block_size[2];
 
302
#ifdef USE_SCALAPACK_DISTR
 
303
    int proc_grid[2];
 
304
#endif
 
305
    double time, gflops;
 
306
    
 
307
    /* create a 2-d GA (global matrix) */
 
308
    dims[0] = matrix_size;
 
309
    dims[1] = matrix_size;
 
310
    block_size[0] = BLOCK_SIZE;
 
311
    block_size[1] = BLOCK_SIZE;
 
312
#ifdef USE_SCALAPACK_DISTR
 
313
    int proc_grid[2];
 
314
    proc_grid[0] = 2;
 
315
    proc_grid[1] = nprocs/2;
 
316
    if(nprocs%2) GA_Error("For ScaLAPACK stle distribution, nprocs must be "
 
317
                         " divisible by 2", 0);
 
318
#endif
 
319
    
 
320
    
 
321
#ifndef BLOCK_CYCLIC
 
322
    g_a = NGA_Create(type, 2, dims, "A", NULL);
 
323
    g_b = GA_Duplicate(g_a, "transposed array B");
 
324
#else
 
325
    g_a = GA_Create_handle();
 
326
    GA_Set_data(g_a, 2, dims, type);
 
327
    GA_Set_array_name(g_a,"A");
 
328
#  ifdef USE_SCALAPACK_DISTR
 
329
    GA_Set_block_cyclic_proc_grid(g_a, block_size, proc_grid);
 
330
#  else
 
331
    GA_Set_block_cyclic(g_a, block_size);    
 
332
#  endif
 
333
    GA_Allocate(g_a);
 
334
    
 
335
    g_b = GA_Create_handle();
 
336
    GA_Set_data(g_b, 2, dims, type);
 
337
    GA_Set_array_name(g_b,"B");
 
338
#  ifdef USE_SCALAPACK_DISTR
 
339
    GA_Set_block_cyclic_proc_grid(g_b, block_size, proc_grid);
 
340
#  else
 
341
    GA_Set_block_cyclic(g_b, block_size);
 
342
#  endif
 
343
    GA_Allocate(g_b);
 
344
    
 
345
#endif
 
346
    
 
347
    /* copy the local matrix into GA */
 
348
    if(me==0) 
 
349
    {
 
350
       lo[0] = 0;
 
351
       hi[0] = matrix_size - 1;
 
352
       lo[1] = 0;
 
353
       hi[1] = matrix_size - 1;
 
354
       ld    = matrix_size;
 
355
       
 
356
       NGA_Put(g_a, lo, hi, A, &ld);
 
357
    }
 
358
    GA_Sync();
 
359
 
 
360
    GA_Transpose(g_a, g_b);
 
361
    time = MP_TIMER();
 
362
    /* The following function does not exist. Not sure what to replace it
 
363
     * with. GA_Lu_solve(char trans, int g_a, int g_b) requiresa an
 
364
     * additioanl GA. */
 
365
    /* GA_Lu('n', g_b); */
 
366
    time = MP_TIMER() - time;
 
367
 
 
368
    /* 2/3 N^3 - 1/2 N^2 flops for LU and 2*N^2 for solver */
 
369
    gflops = ( (((double)matrix_size) * matrix_size)/(time*1.0e+9) *
 
370
               (2.0/3.0 * (double)matrix_size - 0.5) );
 
371
    if(me==0) printf("\nGA_Lu: N=%d flops=%2.5e Gflops, time=%2.5e secs\n\n",
 
372
                     matrix_size, gflops, time);
 
373
 
 
374
#if DEBUG
 
375
    GA_Print(g_a);
 
376
    GA_Print(g_b);
 
377
#endif
 
378
    /* if(me==0) lu(A, matrix_size);     */
 
379
 
 
380
    GA_Destroy(g_a);
 
381
    GA_Destroy(g_b);
 
382
}
 
383
 
 
384
 
 
385
int main(int argc, char **argv) {
 
386
    int i, matrix_size;
 
387
    int heap=20000000, stack=200000000;
 
388
    double *A=NULL;
 
389
    
 
390
    if(argc == 1)
 
391
    {
 
392
      matrix_size = 1234;
 
393
    }
 
394
    else if (argc == 2)
 
395
    {
 
396
      matrix_size = atoi(argv[1]);
 
397
    }
 
398
    else
 
399
    {
 
400
       printf("Usage Error\n\t Usage: <program> [matrix_size]\n");
 
401
       exit(0);
 
402
    }
 
403
 
 
404
    if(matrix_size <= 0) 
 
405
    {
 
406
       printf("Error: matrix size (%d) should be > 0\n",
 
407
              matrix_size);
 
408
       GA_Error("matrix size should be >0", 1);
 
409
    }
 
410
    
 
411
    /* *****************************************************************
 
412
     * Initialize MPI/TCGMSG-MPI, GA and MA
 
413
     * *****************************************************************/
 
414
    MP_INIT(argc,argv);
 
415
    
 
416
    GA_INIT(argc,argv);        /* initialize GA */
 
417
 
 
418
    me     = GA_Nodeid();
 
419
    nprocs = GA_Nnodes();
 
420
 
 
421
    heap /= nprocs;
 
422
    stack /= nprocs;
 
423
 
 
424
    if(! MA_init(MT_F_DBL, stack, heap)) /* initialize MA */
 
425
    {
 
426
       GA_Error("MA_init failed",stack+heap);
 
427
    }
 
428
 
 
429
    /* create/initialize the matrix */    
 
430
    if((A = (double*)malloc(matrix_size*matrix_size*sizeof(double))) == NULL) 
 
431
    {
 
432
       GA_Error("malloc failed", matrix_size*matrix_size*sizeof(double));
 
433
    }
 
434
 
 
435
    for(i=0; i<2; i++) /* 5 runs */
 
436
    {
 
437
       init_array(A, matrix_size);
 
438
#if DEBUG
 
439
       if(me==0) print_array(A, matrix_size);
 
440
#endif
 
441
       /* *****************************************************************
 
442
        * Perform LU Factorization
 
443
        * *****************************************************************/
 
444
       ga_lu(A, matrix_size);
 
445
    }
 
446
    
 
447
    free(A);
 
448
    
 
449
    /* *****************************************************************
 
450
     * Terminate MPI/TCGMSG-MPI, GA and MA
 
451
     * *****************************************************************/
 
452
    if(me==0)printf("Success\n");
 
453
    GA_Terminate();
 
454
 
 
455
    MP_FINALIZE();
 
456
    
 
457
    return 0;
 
458
}
 
459
 
 
460
/**
 
461
 * TODO:
 
462
 *   - LU for non-square matrix
 
463
 *
 
464
 */