← 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/testspd.F

  • 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/testspd.F
1
 
#if HAVE_CONFIG_H
2
 
#   include "config.fh"
3
 
#endif
4
 
c $Id: testspd.F,v 1.8 2004-11-12 22:19:10 edo Exp $
5
 
c***********************************************************************
6
 
c* Source : testllt
7
 
c* Scope  : test LLT SCALAPACK routines
8
 
c*
9
 
c* 04/12/96  GVT  First Implementation
10
 
c*           Giuseppe Vitillaro peppe@unipg.it
11
 
c***********************************************************************
12
 
 
13
 
#define SIZE 100
14
 
#define UPLO 'L'      
15
 
#define TRESH 3.d0
16
 
#define ga_dnormF(g_a) sqrt(ga_ddot(g_a, g_a))
17
 
 
18
 
#define BLOCK_CYCLIC 0
19
 
 
20
 
      
21
 
c***********************************************************************
22
 
c* Program: testllt
23
 
c***********************************************************************
24
 
      program testllt
25
 
      
26
 
      implicit none
27
 
      
28
 
#include "mafdecls.fh"
29
 
#include "global.fh"
30
 
 
31
 
c****
32
 
      integer nproc
33
 
      integer hsize, ssize
34
 
c
35
 
c***  Intitialize a message passing library
36
 
c
37
 
#include "mp3.fh"
38
 
c
39
 
c**** Initialize GA package
40
 
      call ga_initialize()
41
 
c
42
 
c**** get number of nodes and calculate memory allocation      
43
 
      
44
 
      hsize = 6*SIZE*SIZE + 3*SIZE
45
 
      ssize = 3*SIZE*SIZE + 3*SIZE 
46
 
      nproc=ga_nnodes()
47
 
      
48
 
      hsize = (hsize/nproc) + 1
49
 
      ssize = (ssize/nproc) + 1 + 3*256*256
50
 
    
51
 
c**** Initialize MA package
52
 
      if (.not. ma_init(MT_DBL, ssize, hsize))
53
 
     $    call ga_error("ma init failed",ssize+hsize)
54
 
 
55
 
c**** Do Core Tests
56
 
      call ctests(UPLO)
57
 
 
58
 
c**** Exit from the GA package
59
 
      call ga_terminate()
60
 
c
61
 
      call MP_FINALIZE()
62
 
      end
63
 
 
64
 
      
65
 
c***********************************************************************
66
 
c* subroutine: infos
67
 
c*             print informations about this run
68
 
c***********************************************************************      
69
 
      subroutine infos(me, nproc, uplo, eps)
70
 
      implicit none
71
 
#include "mafdecls.fh"
72
 
#include "global.fh"
73
 
c****
74
 
      integer me
75
 
      integer nproc
76
 
      character*1 uplo
77
 
      double precision eps
78
 
c****      
79
 
      
80
 
      if (me.eq.0) then
81
 
         print *, ' '
82
 
         print *, 'Number of nodes  :  ', nproc
83
 
         print *, 'Problem size     :  ', SIZE
84
 
         print *, 'Uplo             :   ', uplo
85
 
         print *, 'Epsilon          :  ', eps
86
 
         call ffflush(6)
87
 
      endif
88
 
 
89
 
      return
90
 
      end
91
 
 
92
 
      
93
 
c***********************************************************************
94
 
c* subroutine: thead
95
 
c*             emit test header output
96
 
c***********************************************************************
97
 
      subroutine thead(header)
98
 
      implicit none
99
 
#include "mafdecls.fh"
100
 
#include "global.fh"
101
 
c****
102
 
      character*(*) header
103
 
c****
104
 
      integer me
105
 
 
106
 
      me = ga_nodeid()
107
 
      if (me.eq.0) then
108
 
         print *, ' '
109
 
         print *, '> ', header, ' '
110
 
         call ffflush(6)
111
 
      endif
112
 
 
113
 
      return
114
 
      end
115
 
 
116
 
      
117
 
c***********************************************************************
118
 
c* subroutine: dthtest
119
 
c*             test a double precision against the THRESH parameter
120
 
c***********************************************************************
121
 
      subroutine dthtest(msg, dval)
122
 
      implicit none
123
 
#include "mafdecls.fh"
124
 
#include "global.fh"
125
 
c****      
126
 
      character*(*) msg
127
 
      double precision dval
128
 
c****
129
 
      integer me
130
 
 
131
 
      me = ga_nodeid()
132
 
 
133
 
      if (me.eq.0) then
134
 
         print *, '  ', msg, dval, ' '
135
 
         if (dval.lt.TRESH) then
136
 
            print *, '> success '
137
 
         else
138
 
            print *, '> failure '
139
 
         endif
140
 
         call ffflush(6)
141
 
      endif
142
 
      
143
 
      return
144
 
      end
145
 
 
146
 
      
147
 
c***********************************************************************
148
 
c* subrotine: stest
149
 
c*            test solver result
150
 
c***********************************************************************
151
 
      subroutine stest(irc, ierc)
152
 
      implicit none
153
 
#include "mafdecls.fh"
154
 
#include "global.fh"
155
 
c****      
156
 
      integer irc
157
 
      integer ierc
158
 
c****
159
 
      integer me
160
 
 
161
 
      me = ga_nodeid()
162
 
 
163
 
      if (me.eq.0) then
164
 
         
165
 
         if (irc.eq.0) then
166
 
            print *, '> LLT '
167
 
         elseif (irc.gt.0) then
168
 
            print *, '> LU  '
169
 
         endif
170
 
 
171
 
         if (irc.eq.ierc .or. (irc.gt.0 .and. ierc.gt.0)) then
172
 
            print *, '> success: expected factorization'
173
 
         else
174
 
            print *, '> failure: not expected factorization',irc,ierc
175
 
         endif
176
 
         
177
 
         call ffflush(6)
178
 
         
179
 
      endif
180
 
      
181
 
      return
182
 
      end
183
 
 
184
 
      
185
 
c***********************************************************************
186
 
c* subroutine: emsg
187
 
c*             error message
188
 
c***********************************************************************
189
 
      subroutine emsg(msg, ival)
190
 
      implicit none
191
 
#include "mafdecls.fh"
192
 
#include "global.fh"
193
 
c****      
194
 
      character*(*) msg
195
 
      integer ival
196
 
c****
197
 
      integer me
198
 
 
199
 
      me = ga_nodeid()
200
 
 
201
 
      if (me.eq.0) then
202
 
         print *, '  >>> ', msg, ival, ' '
203
 
         print *, '> failure '
204
 
         call ffflush(6)
205
 
      endif
206
 
      
207
 
      return
208
 
      end
209
 
 
210
 
      
211
 
 
212
 
      
213
 
c***********************************************************************
214
 
c* subroutine: ctests
215
 
c*             do coretests for LLT Cholesky factorization, solver
216
 
c***********************************************************************
217
 
      subroutine ctests(uplo)
218
 
      
219
 
      implicit none
220
 
            
221
 
#include "mafdecls.fh"
222
 
#include "global.fh"
223
 
      integer ga_cholesky
224
 
      external ga_cholesky
225
 
c****
226
 
      character*1 uplo
227
 
c****      
228
 
      integer n
229
 
      parameter (n=SIZE)
230
 
 
231
 
      external ga_llt_f, ga_llt_s, ga_llt_i
232
 
      integer ga_llt_f, ga_llt_i
233
 
 
234
 
      double precision A(n,n), X(n,n)
235
 
      
236
 
      integer g_A, g_B, g_X, g_Y
237
 
      integer g_A1, g_X1, g_Y1, g_Y2
238
 
#if BLOCK_CYCLIC
239
 
      integer g_AA, g_BB, g_CC, g_YY, g_ZZ
240
 
      integer dims(2), proc_grid(2), block(2)
241
 
      integer g1, g2
242
 
#endif
243
 
 
244
 
      integer hsA
245
 
 
246
 
      integer i, j
247
 
      integer nproc, me
248
 
      double precision buf(n)
249
 
      
250
 
      double precision dnA, dnX, dnX1, dnY1, dnS, dnF, dnI
251
 
      integer irc
252
 
      
253
 
      double precision eps,  dlamch
254
 
 
255
 
      eps   =  dlamch('p')
256
 
      
257
 
      nproc = ga_nnodes()
258
 
      me    = ga_nodeid()
259
 
 
260
 
      call infos(me, nproc, uplo, eps)
261
 
 
262
 
      
263
 
c****************************************
264
 
c* Initialize tests variables
265
 
c* Generate a Lower triangula matrix
266
 
c* with large positive diagonal elements
267
 
c* so the LU decomposition will be
268
 
c* just this matrix      
269
 
c****************************************
270
 
c**** Initialize local arrays A, X
271
 
c**** they are local copies of the corrispondent
272
 
c**** global arrays
273
 
      do j = 1, n
274
 
         do i = 1, n
275
 
            X(i,j) = dsin(1.d0 * (i * j))
276
 
            A(i,j) = 0.d0
277
 
            if (i.gt.j) then
278
 
               A(i,j) = dsin(1.d0 * (i + j))
279
 
            endif
280
 
            if (i.eq.j) then
281
 
               A(i,j) = SIZE*dabs(dsin(10.d0 * i))
282
 
            endif
283
 
         end do ! i
284
 
      end do    ! j
285
 
      
286
 
c**** Create global arrays
287
 
#if BLOCK_CYCLIC
288
 
      if (me.eq.0) then
289
 
        write(6,*) '*'
290
 
        write(6,*) '* Creating Block-Cyclic Arrays'
291
 
        write(6,*) '*'
292
 
      endif
293
 
      dims(1) = n
294
 
      dims(2) = n
295
 
      block(1) = 64
296
 
      block(2) = 64
297
 
      call factor(nproc,g1,g2)
298
 
      proc_grid(1) = g1
299
 
      proc_grid(2) = g2
300
 
 
301
 
      g_A = ga_create_handle()
302
 
      call ga_set_data(g_A,2,dims,MT_DBL)
303
 
      call ga_set_block_cyclic_proc_grid(g_A,block,proc_grid)
304
 
      if (.not.ga_allocate(g_A)) 
305
 
     &     call ga_error(' ga_create A failed ', 2)
306
 
 
307
 
      g_B = ga_create_handle()
308
 
      call ga_set_data(g_B,2,dims,MT_DBL)
309
 
      call ga_set_block_cyclic_proc_grid(g_B,block,proc_grid)
310
 
      if (.not.ga_allocate(g_B)) 
311
 
     &     call ga_error(' ga_create B failed ', 2)
312
 
 
313
 
      g_A1 = ga_create_handle()
314
 
      call ga_set_data(g_A1,2,dims,MT_DBL)
315
 
      call ga_set_block_cyclic_proc_grid(g_A1,block,proc_grid)
316
 
      if (.not.ga_allocate(g_A1)) 
317
 
     &     call ga_error(' ga_create A1 failed ', 2)
318
 
 
319
 
      g_X = ga_create_handle()
320
 
      call ga_set_data(g_X,2,dims,MT_DBL)
321
 
      call ga_set_block_cyclic_proc_grid(g_X,block,proc_grid)
322
 
      if (.not.ga_allocate(g_X)) 
323
 
     &     call ga_error(' ga_create X failed ', 2)
324
 
 
325
 
      g_X1 = ga_create_handle()
326
 
      call ga_set_data(g_X1,2,dims,MT_DBL)
327
 
      call ga_set_block_cyclic_proc_grid(g_X1,block,proc_grid)
328
 
      if (.not.ga_allocate(g_X1)) 
329
 
     &     call ga_error(' ga_create X1 failed ', 2)
330
 
 
331
 
      g_AA = ga_create_handle()
332
 
      call ga_set_data(g_AA,2,dims,MT_DBL)
333
 
      if (.not.ga_allocate(g_AA)) 
334
 
     &     call ga_error(' ga_create AA failed ', 2)
335
 
 
336
 
      g_BB = ga_create_handle()
337
 
      call ga_set_data(g_BB,2,dims,MT_DBL)
338
 
      if (.not.ga_allocate(g_BB)) 
339
 
     &     call ga_error(' ga_create BB failed ', 2)
340
 
 
341
 
      g_CC = ga_create_handle()
342
 
      call ga_set_data(g_CC,2,dims,MT_DBL)
343
 
      if (.not.ga_allocate(g_CC)) 
344
 
     &     call ga_error(' ga_create CC failed ', 2)
345
 
 
346
 
      dims(2) = 1
347
 
 
348
 
      g_Y = ga_create_handle()
349
 
      call ga_set_data(g_Y,2,dims,MT_DBL)
350
 
      call ga_set_block_cyclic_proc_grid(g_Y,block,proc_grid)
351
 
      if (.not.ga_allocate(g_Y)) 
352
 
     &     call ga_error(' ga_create Y failed ', 2)
353
 
 
354
 
      g_Y1 = ga_create_handle()
355
 
      call ga_set_data(g_Y1,2,dims,MT_DBL)
356
 
      call ga_set_block_cyclic_proc_grid(g_Y1,block,proc_grid)
357
 
      if (.not.ga_allocate(g_Y1)) 
358
 
     &     call ga_error(' ga_create Y1 failed ', 2)
359
 
 
360
 
      g_Y2 = ga_create_handle()
361
 
      call ga_set_data(g_Y2,2,dims,MT_DBL)
362
 
      call ga_set_block_cyclic_proc_grid(g_Y2,block,proc_grid)
363
 
      if (.not.ga_allocate(g_Y2)) 
364
 
     &     call ga_error(' ga_create Y1 failed ', 2)
365
 
 
366
 
      g_YY = ga_create_handle()
367
 
      call ga_set_data(g_YY,2,dims,MT_DBL)
368
 
      if (.not.ga_allocate(g_YY)) 
369
 
     &     call ga_error(' ga_create YY failed ', 2)
370
 
 
371
 
      g_ZZ = ga_create_handle()
372
 
      call ga_set_data(g_ZZ,2,dims,MT_DBL)
373
 
      if (.not.ga_allocate(g_ZZ)) 
374
 
     &     call ga_error(' ga_create YY failed ', 2)
375
 
#else
376
 
      if (.not. ga_create(MT_DBL, n, n, 'A', 1, 1, g_A))
377
 
     &     call ga_error(' ga_create A failed ', 2)
378
 
      
379
 
      if (.not. ga_create(MT_DBL, n, n, 'B', 1, 1, g_B))
380
 
     &     call ga_error(' ga_create B failed ', 2)
381
 
 
382
 
      if (.not. ga_create(MT_DBL, n, n, 'A1', 1, 1, g_A1))
383
 
     &     call ga_error(' ga_create A1 failed ', 2)
384
 
 
385
 
      if (.not. ga_create(MT_DBL, n, n, 'X', 1, 1, g_X))
386
 
     &     call ga_error(' ga_create X failed ', 2)
387
 
 
388
 
      if (.not. ga_create(MT_DBL, n, n, 'X1', 1, 1, g_X1))
389
 
     &     call ga_error(' ga_create X1 failed ', 2)
390
 
      
391
 
      if (.not. ga_create(MT_DBL, n, 1, 'Y', 1, 1, g_Y))
392
 
     &     call ga_error(' ga_create Y failed ', 2)
393
 
      
394
 
      if (.not. ga_create(MT_DBL, n, 1, 'Y1', 1, 1, g_Y1))
395
 
     &     call ga_error(' ga_create Y1 failed ', 2)
396
 
 
397
 
      if (.not. ga_create(MT_DBL, n, 1, 'Y2', 1, 1, g_Y2))
398
 
     &     call ga_error(' ga_create Y2 failed ', 2)
399
 
#endif
400
 
 
401
 
c**** Fill in arrays A, X
402
 
      do j = me+1, n, nproc
403
 
         call ga_put(g_A, 1, n, j, j, A(1,j), n)
404
 
         call ga_put(g_X, 1, n, j, j, X(1,j), n)
405
 
      end do    ! j 
406
 
 
407
 
c**** Create A, B, Y
408
 
c**** It is granted that A, B, Y will not change
409
 
c**** during the execution of the test LLT driver      
410
 
      
411
 
c**** B = A * A^
412
 
#if BLOCK_CYCLIC
413
 
      call ga_copy(g_A, g_AA)
414
 
      call ga_dgemm('N', 'T', n, n, n, 1.d0, g_AA, g_AA, 0.d0, g_BB)
415
 
      call ga_copy(g_BB, g_B)
416
 
#else
417
 
      call ga_dgemm('N', 'T', n, n, n, 1.d0, g_A, g_A, 0.d0, g_B)
418
 
#endif
419
 
 
420
 
c**** A = B = A * A^
421
 
      call ga_copy(g_B, g_A)
422
 
      
423
 
c**** B = A * X      
424
 
#if BLOCK_CYCLIC
425
 
      call ga_copy(g_A, g_AA)
426
 
      call ga_copy(g_X, g_CC)
427
 
      call ga_dgemm('N', 'N', n, n, n, 1.d0, g_AA, g_CC, 0.d0, g_BB)
428
 
      call ga_copy(g_BB,g_B)
429
 
#else
430
 
      call ga_dgemm('N', 'N', n, n, n, 1.d0, g_A, g_X, 0.d0, g_B)
431
 
#endif
432
 
 
433
 
c**** Copy B(1:n,1:1) to Y(1:n) so Y will X column 1
434
 
      call ga_copy_patch('N', g_B, 1, n, 1, 1, g_Y, 1, n, 1, 1)
435
 
 
436
 
c**** dnA = ||A||
437
 
      dnA = ga_dnormF(g_A)
438
 
 
439
 
      
440
 
c****************************************************
441
 
c Test Cholesky factorization external interface
442
 
c****************************************************
443
 
      call thead("Test Cholesky factorization")
444
 
 
445
 
c**** copy A in X
446
 
      call ga_copy(g_A, g_X)
447
 
 
448
 
c**** call Cholesky factorization routine in ScaLAPACK PDPOTRF
449
 
c**** to obtain an LL'/U'U factorization (external interface)
450
 
      irc = ga_cholesky(uplo, g_X)
451
 
      
452
 
c**** if return code 'zero' is OK
453
 
      if (irc.eq.0) then
454
 
 
455
 
c****    A1 = X * X^ or A1 = X^ * X
456
 
         if (uplo.eq.'L') then
457
 
#if BLOCK_CYCLIC
458
 
          call ga_copy(g_X, g_CC)
459
 
          call ga_dgemm('N', 'T', n, n, n, 1.d0, g_CC, g_CC, 0.d0, g_AA)
460
 
          call ga_copy(g_AA, g_A1)
461
 
#else
462
 
          call ga_dgemm('N', 'T', n, n, n, 1.d0, g_X, g_X, 0.d0, g_A1)
463
 
#endif
464
 
         else
465
 
#if BLOCK_CYCLIC
466
 
          call ga_copy(g_X, g_CC)
467
 
          call ga_dgemm('T', 'N', n, n, n, 1.d0, g_CC, g_CC, 0.d0, g_AA)
468
 
          call ga_copy(g_AA, g_A1)
469
 
#else
470
 
          call ga_dgemm('T', 'N', n, n, n, 1.d0, g_X, g_X, 0.d0, g_A1)
471
 
#endif
472
 
         endif
473
 
         
474
 
c****    A1 = A - A1 = A - X * X^ or A1 = A - A1 = A - X^ * X
475
 
         call ga_add(1.0d0, g_A1, -1.0d0, g_A, g_A1)
476
 
      
477
 
c****    dnF = ||A - A1*A1^|| / (||A|| * N * eps) : SHAPE L
478
 
c****    dnF = ||A - A1^*A1|| / (||A|| * N * eps) : SHAPE U         
479
 
         dnF = ga_dnormF(g_A1) / (dnA * n * eps)
480
 
 
481
 
         if (uplo.eq.'L') then
482
 
            call dthtest("||LL' - A|| / (||A|| * N * eps) =", dnF)
483
 
         else   
484
 
            call dthtest("||U'U - A|| / (||A|| * N * eps) =", dnF)
485
 
         endif
486
 
         
487
 
c**** if return code is > 0
488
 
      else
489
 
         call emsg('It is not positive definite the minor of order:',
490
 
     &             irc)
491
 
      endif
492
 
 
493
 
c****************************************************
494
 
c Test Cholesky factorization and solver
495
 
c internal interfaces with a NxN RHS array
496
 
c****************************************************
497
 
      call thead("Test II Cholesky solver with a NxN RHS array")
498
 
 
499
 
c**** call Cholesky factorization routine in ScaLAPACK PDPOTRF
500
 
c**** to obtain an LL'/U'U factorization
501
 
c**** (internal interface: it will not destroy A)
502
 
#if BLOCK_CYCLIC
503
 
         call ga_copy(g_A, g_AA)
504
 
#endif
505
 
      hsA = 0
506
 
      irc = ga_llt_f(uplo, g_A, hsA)
507
 
      
508
 
c**** if return code from ga_llt_f is zero
509
 
      if (irc.eq.0) then
510
 
         
511
 
c****    copy B in X
512
 
         call ga_copy(g_B, g_X)
513
 
         
514
 
c****    call Cholesky solver routine in ScaLAPACK PDPOTRS
515
 
c****    internal interface with an NxN RHS GA
516
 
         call ga_llt_s(uplo, g_A, g_X, hsA)
517
 
 
518
 
c****    A1 = A * X
519
 
#if BLOCK_CYCLIC
520
 
         call ga_copy(g_X, g_BB)
521
 
         call ga_dgemm('N', 'N', n, n, n, 1.d0,
522
 
     &        g_AA, g_BB, 0.d0, g_CC)
523
 
         call ga_copy(g_CC,g_A1)
524
 
#else
525
 
         call ga_dgemm('N', 'N', n, n, n, 1.d0,
526
 
     &        g_A, g_X, 0.d0, g_A1)
527
 
#endif
528
 
 
529
 
c****    A1 = A1 - B = A * X - B
530
 
         call ga_add(1.d0, g_A1, -1.d0, g_B, g_A1)
531
 
         
532
 
c****    dnX = ||X||
533
 
         dnX = ga_dnormF(g_X)
534
 
         
535
 
c****    dnS = ||AX - B|| / (||A|| * ||X|| * N * eps)
536
 
         dnS = ga_dnormF(G_A1) / (dnA * dnX * n * eps)
537
 
 
538
 
         call dthtest("||A*X - B||/(||A||*||X||*n*eps) =", dnS)
539
 
 
540
 
c**** if return code is > 0
541
 
      else
542
 
         call emsg('It is not positive definite the minor of order:',
543
 
     &        irc)
544
 
      endif
545
 
 
546
 
 
547
 
c****************************************************
548
 
c Test Cholesky factorization and solver
549
 
c internal interfaces with a single vector RHS
550
 
c****************************************************
551
 
      call thead("Test II Cholesky solver with a single vector RHS")
552
 
 
553
 
c**** call Cholesky factorization routine in ScaLAPACK PDPOTRF
554
 
c**** to obtain an LL'/U'U factorization
555
 
c**** (internal interface: it will not destroy A)
556
 
#if BLOCK_CYCLIC
557
 
         call ga_copy(g_AA, g_A)
558
 
#endif
559
 
      hsA = 0
560
 
      irc = ga_llt_f(uplo, g_A, hsA)
561
 
      
562
 
c**** if return code from ga_llt_f is zero
563
 
      if (irc.eq.0) then
564
 
         
565
 
c****    copy Y in Y1
566
 
         call ga_copy(g_Y, g_Y1)
567
 
         
568
 
c****    call Cholesky solver routine in ScaLAPACK PDPOTRS
569
 
c****    internal interface with a single vector RHS
570
 
         call ga_llt_s(uplo, g_A, g_Y1, hsA)
571
 
 
572
 
c****    dnY1 = ||Y1||
573
 
         dnY1 = ga_dnormF(g_Y1)
574
 
                  
575
 
c****    Y1 = A * Y1
576
 
#if BLOCK_CYCLIC
577
 
         call ga_copy(g_Y1, g_YY)
578
 
         call ga_dgemm('N', 'N', n, 1, n, 1.d0,
579
 
     &        g_AA, g_YY, 0.d0, g_ZZ)
580
 
         call ga_copy(g_ZZ,g_Y2)
581
 
#else
582
 
         call ga_dgemm('N', 'N', n, 1, n, 1.d0,
583
 
     &        g_A, g_Y1, 0.d0, g_Y2)
584
 
#endif
585
 
 
586
 
c****    Y1 = Y1 - Y = A * Y1 - Y
587
 
         call ga_add(1.d0, g_Y2, -1.d0, g_Y, g_Y1)
588
 
 
589
 
c****    dnS = ||AY1 - Y|| / (||A|| * ||Y1|| * N * eps)
590
 
         dnS = ga_dnormF(G_A1) / (dnA * dnY1 * n * eps)
591
 
 
592
 
         call dthtest("||A*X - V||/(||A||*||X||*n*eps) =", dnS)
593
 
 
594
 
c**** if return code is > 0
595
 
      else
596
 
         call emsg('It is not positive definite the minor of order:',
597
 
     &        irc)         
598
 
      endif
599
 
      
600
 
c****************************************************
601
 
c Test Cholesky solver with a NxN RHS array
602
 
c****************************************************
603
 
      call thead("Test EI Cholesky solver with a NxN RHS array")
604
 
 
605
 
c**** copy B in X
606
 
      call ga_copy(g_B, g_X)
607
 
         
608
 
c**** call Cholesky L/U solver with a NxN RHS
609
 
#if BLOCK_CYCLIC
610
 
      call ga_copy(g_AA,g_A)
611
 
#endif
612
 
      irc = ga_llt_solve(g_A, g_X)
613
 
      
614
 
c**** if return code from ga_llt_solve is zero
615
 
      if (irc.eq.0) then
616
 
         
617
 
c****    A1 = A * X
618
 
#if BLOCK_CYCLIC
619
 
         call ga_copy(g_X,g_BB)
620
 
         call ga_dgemm('N', 'N', n, n, n, 1.d0,
621
 
     &        g_AA, g_BB, 0.d0, g_CC)
622
 
         call ga_copy(g_CC,g_A1)
623
 
#else
624
 
         call ga_dgemm('N', 'N', n, n, n, 1.d0,
625
 
     &        g_A, g_X, 0.d0, g_A1)
626
 
#endif
627
 
 
628
 
c****    A1 = A1 - B = A * X - B
629
 
         call ga_add(1.d0, g_A1, -1.d0, g_B, g_A1)
630
 
 
631
 
c****    dnX = ||X||
632
 
         dnX = ga_dnormF(g_X)
633
 
         
634
 
c****    dnS = ||AX - B|| / (||A|| * ||X|| * N * eps)
635
 
         dnS = ga_dnormF(G_A1) / (dnA * dnX * n * eps)
636
 
 
637
 
         call dthtest("||A*X - B||/(||A||*||X||*n*eps) =", dnS)
638
 
         
639
 
c**** if return code is > 0
640
 
      else
641
 
         call emsg('It is not positive definite the minor of order:',
642
 
     &        irc)                  
643
 
      endif
644
 
 
645
 
      
646
 
c****************************************************
647
 
c Test Cholesky solver with a single vector RHS
648
 
c****************************************************
649
 
      call thead("Test EI Cholesky solver with a single vector RHS")
650
 
      
651
 
c**** copy Y in Y1
652
 
      call ga_copy(g_Y, g_Y1)
653
 
 
654
 
c**** call Cholesky L/U solver with a single vector RHS
655
 
#if BLOCK_CYCLIC
656
 
      call ga_copy(g_AA,g_A)
657
 
#endif
658
 
      irc = ga_llt_solve(g_A, g_Y1)
659
 
      
660
 
c**** if return code from ga_llt_solve is zero
661
 
      if (irc.eq.0) then
662
 
 
663
 
c****    dnY1 = ||Y1||
664
 
         dnY1 = ga_dnormF(g_Y1)
665
 
                  
666
 
c****    Y1 = A * Y1
667
 
#if BLOCK_CYCLIC
668
 
         call ga_copy(g_Y1,g_YY)
669
 
         call ga_dgemm('N', 'N', n, 1, n, 1.d0,
670
 
     &        g_AA, g_YY, 0.d0, g_ZZ)
671
 
         call ga_copy(g_ZZ,g_Y2)
672
 
#else
673
 
         call ga_dgemm('N', 'N', n, 1, n, 1.d0,
674
 
     &        g_A, g_Y1, 0.d0, g_Y2)
675
 
#endif
676
 
 
677
 
c****    Y1 = Y1 - Y = A * Y1 - Y
678
 
         call ga_add(1.d0, g_Y2, -1.d0, g_Y, g_Y1)
679
 
 
680
 
c****    dnS = ||AY1 - Y|| / (||A|| * ||Y1|| * N * eps)
681
 
         dnS = ga_dnormF(G_A1) / (dnA * dnY1 * n * eps)
682
 
         
683
 
         call dthtest("||A*X - V||/(||A||*||X||*n*eps) =", dnS)
684
 
         
685
 
c**** if return code is > 0
686
 
      else
687
 
         call emsg('It is not positive definite the minor of order:',
688
 
     &        irc)                  
689
 
      endif
690
 
 
691
 
      
692
 
c****************************************************
693
 
c Test Cholesky factorization and inversion
694
 
c internal interfaces 
695
 
c****************************************************
696
 
      call thead("Test II inversion of an SPD matrix")
697
 
 
698
 
c**** copy A in X
699
 
#if BLOCK_CYCLIC
700
 
      call ga_copy(g_AA, g_A)
701
 
#endif
702
 
      call ga_copy(g_A, g_X)
703
 
 
704
 
c**** call Cholesky factorization routine in ScaLAPACK PDPOTRF
705
 
c**** to obtain an LL'/U'U factorization
706
 
c**** (internal interface: it will not destroy X1)
707
 
      hsA = 0
708
 
      irc = ga_llt_f(uplo, g_X, hsA)
709
 
      
710
 
c**** if return code from ga_llt_f is zero
711
 
      if (irc.eq.0) then
712
 
         
713
 
c****    call Cholesky inversion routine in ScaLAPACK PDPOTRI
714
 
c****    internal interface
715
 
         irc = ga_llt_i(uplo, g_X, hsA)
716
 
         
717
 
c****    if the inversion could be done         
718
 
         if (irc.eq.0) then
719
 
 
720
 
c****       dnX = ||X|| = ||invA||
721
 
            dnX = ga_dnormF(g_X)
722
 
 
723
 
c****       A1 = A * X = A * invA
724
 
#if BLOCK_CYCLIC
725
 
            call ga_copy(g_X,g_BB)
726
 
            call ga_dgemm('N', 'N', n, n, n, 1.d0,
727
 
     &           g_AA, g_BB, 0.d0, g_CC)
728
 
            call ga_copy(g_CC,g_A1)
729
 
#else
730
 
            call ga_dgemm('N', 'N', n, n, n, 1.d0,
731
 
     &           g_A, g_X, 0.d0, g_A1)
732
 
#endif
733
 
 
734
 
c****       A1 = A1 - I = A * invA - I
735
 
            do j = me+1, n, nproc
736
 
               call ga_get(g_A1, j, j, j, j, buf, 1)
737
 
               buf(1) = buf(1) - 1.d0
738
 
               call ga_put(g_A1, j, j, j, j, buf, 1)
739
 
            end do ! j
740
 
 
741
 
c****       dnI = ||A*invA - I|| / (||A| * ||invA|| * N * eps)
742
 
            dnI = ga_dnormF(g_A1) / (dnA * dnX * n * eps)
743
 
 
744
 
            call dthtest('||A*invA - I||/(||A||*||invA||*n*eps) =',
745
 
     &           dnI)
746
 
         
747
 
         else
748
 
c****       otherwise if the ga_llt_i return code is > 0
749
 
            call emsg('there is a zero diagonal element:', irc)
750
 
         endif
751
 
      
752
 
c**** if return code is > 0
753
 
      else
754
 
         call emsg('It is not positive definite the minor of order:',
755
 
     &        irc)                  
756
 
      endif
757
 
 
758
 
      
759
 
c****************************************************
760
 
c Test Cholesky inversion
761
 
c****************************************************
762
 
      call thead("Test inversion of an SPD matrix")
763
 
 
764
 
c**** copy A in X
765
 
#if BLOCK_CYCLIC
766
 
      call ga_copy(g_AA, g_A)
767
 
#endif
768
 
      call ga_copy(g_A, g_X)
769
 
 
770
 
c**** call Cholesky L/U inversion
771
 
      irc = ga_spd_invert(g_X)
772
 
      
773
 
c**** if return code from ga_spd_invert is zero
774
 
      if (irc.eq.0) then
775
 
 
776
 
c****    dnX = ||X|| = ||invA||
777
 
         dnX = ga_dnormF(g_X)
778
 
         
779
 
c****    A1 = A * X = A * invA
780
 
#if BLOCK_CYCLIC
781
 
         call ga_copy(g_X, g_BB)
782
 
         call ga_dgemm('N', 'N', n, n, n, 1.d0,
783
 
     &        g_AA, g_BB, 0.d0, g_CC)
784
 
         call ga_copy(g_CC, g_A1)
785
 
#else
786
 
         call ga_dgemm('N', 'N', n, n, n, 1.d0,
787
 
     &        g_A, g_X, 0.d0, g_A1)
788
 
#endif
789
 
 
790
 
c****    A1 = A1 - I = A * invA - I
791
 
         do j = me+1, n, nproc
792
 
            call ga_get(g_A1, j, j, j, j, buf, 1)
793
 
            buf(1) = buf(1) - 1.d0
794
 
            call ga_put(g_A1, j, j, j, j, buf, 1)
795
 
         end do ! j
796
 
 
797
 
c****    dnI = ||A*invA - I|| / (||A|| * ||invA|| * N * eps)
798
 
         dnI = ga_dnormF(g_A1) / (dnA * dnX * n * eps)
799
 
 
800
 
         call dthtest('||A*invA - I||/(||A||*||invA||*n*eps) =',
801
 
     &        dnI)
802
 
         
803
 
c**** if return code is < 0
804
 
      elseif (irc.lt.0) then
805
 
         call emsg('there is a zero diagonal element:', irc)
806
 
c**** if return code is > 0         
807
 
      elseif (irc.gt.0) then
808
 
         call emsg('It is not positive definite the minor of order:',
809
 
     &        irc)
810
 
      endif
811
 
 
812
 
      
813
 
c****************************************************
814
 
c* Test general solver with a NxN RSH
815
 
c* simmetric positive definite array       
816
 
c****************************************************
817
 
      call thead(
818
 
     & 'Test solver for a symmetric P.D. matrix and NxN RHS'
819
 
     &     )
820
 
      
821
 
c**** copy B in X
822
 
      call ga_copy(g_B, g_X)
823
 
         
824
 
c**** call general solver with a NxN RHS p.d. symmetric array
825
 
#if BLOCK_CYCLIC
826
 
      call ga_copy(g_AA, g_A)
827
 
#endif
828
 
      irc = ga_solve(g_A, g_X)
829
 
 
830
 
c**** A1 = A * X
831
 
#if BLOCK_CYCLIC
832
 
      call ga_copy(g_X, g_BB)
833
 
      call ga_dgemm('N', 'N', n, n, n, 1.d0,
834
 
     &     g_AA, g_BB, 0.d0, g_CC)
835
 
      call ga_copy(g_CC, g_A1)
836
 
#else
837
 
      call ga_dgemm('N', 'N', n, n, n, 1.d0,
838
 
     &     g_A, g_X, 0.d0, g_A1)
839
 
#endif
840
 
 
841
 
c**** A1 = A1 - B = A * X - B
842
 
      call ga_add(1.d0, g_A1, -1.d0, g_B, g_A1)
843
 
 
844
 
c**** dnX = ||X||
845
 
      dnX = ga_dnormF(g_X)
846
 
      
847
 
c**** dnS = ||AX - B|| / (||A|| * ||X|| * N * eps)
848
 
      dnS = ga_dnormF(G_A1) / (dnA * dnX * n * eps)
849
 
 
850
 
      call dthtest("||A*X - B||/(||A||*||X||*n*eps) =", dnS)
851
 
 
852
 
      call stest(irc, 0)
853
 
 
854
 
      
855
 
c****************************************************
856
 
c* Test general solver with a NxN RSH
857
 
c* simmetric non positive definite array       
858
 
c****************************************************
859
 
      call thead(
860
 
     & 'Test solver for a symmetric NON P.D. matrix and NxN RHS')
861
 
      
862
 
c**** copy B on X1
863
 
      call ga_copy(g_B, g_X1)
864
 
      
865
 
c**** and now symmetrize X1
866
 
c**** so we obtain a symmetric matrix in X1
867
 
      call ga_symmetrize(g_X1)
868
 
 
869
 
c**** copy B in X
870
 
      call ga_copy(g_B, g_X)
871
 
      
872
 
c**** call general solver with a NxN RHS symmetric array
873
 
#if BLOCK_CYCLIC
874
 
      call ga_copy(g_X1, g_BB)
875
 
      irc = ga_solve(g_BB, g_X)
876
 
#else
877
 
      irc = ga_solve(g_X1, g_X)
878
 
#endif
879
 
      
880
 
c**** dnX1 = ||X1||
881
 
      dnX1 = ga_dnormF(g_X1)
882
 
 
883
 
c**** dnX = ||X||
884
 
      dnX = ga_dnormF(g_X)
885
 
 
886
 
c**** A1 = X1 * X
887
 
#if BLOCK_CYCLIC
888
 
      call ga_copy(g_X1,g_AA)
889
 
      call ga_copy(g_X,g_BB)
890
 
      call ga_dgemm('N', 'N', n, n, n, 1.d0,
891
 
     &     g_AA, g_BB, 0.d0, g_CC)
892
 
      call ga_copy(g_CC,g_A1)
893
 
#else
894
 
      call ga_dgemm('N', 'N', n, n, n, 1.d0,
895
 
     &     g_X1, g_X, 0.d0, g_A1)
896
 
#endif
897
 
 
898
 
c**** A1 = A1 - B = X1 * X - B
899
 
      call ga_add(1.d0, g_A1, -1.d0, g_B, g_A1)
900
 
 
901
 
c**** dnS = ||X1*X - B|| / (||X1|| * ||X|| * N * eps)
902
 
      dnS = ga_dnormF(G_A1) / (dnX1 * dnX * n * eps)
903
 
 
904
 
      call dthtest("||A*X - B||/(||A||*||X||*n*eps) =", dnS)
905
 
     
906
 
      call stest(irc, 1)
907
 
 
908
 
      
909
 
c****************************************************
910
 
c CTESTS exit code
911
 
c****************************************************      
912
 
c**** just print a newline and return      
913
 
      if (me.eq.0) then
914
 
         print *, ' '
915
 
      endif
916
 
 
917
 
      return
918
 
 
919
 
      end
920
 
#if BLOCK_CYCLIC
921
 
      subroutine factor(p,idx,idy)
922
 
      implicit none
923
 
      integer i,j,p,idx,idy,it
924
 
      integer ip,ifac,pmax,prime(1280)
925
 
      integer fac(1280)
926
 
c
927
 
      i = 1
928
 
      ip = p
929
 
c
930
 
c    factor p completely
931
 
c    first, find all prime numbers less than or equal to p
932
 
c
933
 
      pmax = 0
934
 
      do i = 2, p
935
 
        do j = 1, pmax
936
 
          if (mod(i,prime(j)).eq.0) go to 100
937
 
        end do
938
 
        pmax = pmax + 1
939
 
        prime(pmax) = i
940
 
  100   continue
941
 
      end do
942
 
c
943
 
c    find all prime factors of p
944
 
c
945
 
      ifac = 0
946
 
      do i = 1, pmax
947
 
  200   if (mod(ip,prime(i)).eq.0) then
948
 
          ifac = ifac + 1
949
 
          fac(ifac) = prime(i)
950
 
          ip = ip/prime(i)
951
 
          go to 200
952
 
        endif
953
 
      end do
954
 
c
955
 
c    determine two factors of p of approximately the
956
 
c    same size
957
 
c
958
 
      idx = 1
959
 
      idy = 1
960
 
      do i = ifac, 1, -1
961
 
        if (idx.le.idy) then
962
 
          idx = fac(i)*idx
963
 
        else
964
 
          idy = fac(i)*idy
965
 
        endif
966
 
      end do
967
 
      return
968
 
      end
969
 
#endif