← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/tools/ga-4-3/global/testing/testspd.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

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