← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/util/ga_lkain_2cpl3.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/util/ga_lkain_2cpl3.F
1
 
      subroutine ga_lkain_2cpl3(rtdb, g_x, g_b, g_x_im, g_b_im,
2
 
     &   product, precond, 
3
 
     $   tol, mmaxsub, maxiter, odiff, oprint, omega, limag,
4
 
     &   lifetime, gamwidth, ncomp)
 
1
      subroutine ga_lkain_2cpl3(rtdb, 
 
2
     &                          g_x, g_b, 
 
3
     &                          g_x_im, g_b_im,
 
4
     &                          product, precond, 
 
5
     $                          tol, mmaxsub, maxiter, 
 
6
     &                          odiff, oprint, omega, limag,
 
7
     &                          lifetime, gamwidth, ncomp)
5
8
 
6
 
c     $Id: ga_lkain_2cpl3.F 19707 2010-10-29 17:59:36Z d3y133 $
 
9
c     $Id: ga_lkain_2cpl3.F 23357 2013-01-04 18:00:34Z niri $
7
10
 
8
11
      implicit none
9
12
#include "errquit.fh"
30
33
      integer maxiter           ! [input] maximum no. of iterations
31
34
      logical odiff             ! [input] use differences in product
32
35
      logical oprint            ! [input] print flag
 
36
 
 
37
      integer ipm
33
38
c
34
39
c     Solves the linear equations A(X)=0 for multiple vectors.
35
40
c
61
66
c
62
67
c     Needs to be extended to store the sub-space vectors out-of-core
63
68
c     at least while the product() routine is being executed.
64
 
c
65
 
      integer iter, n, n2, nvec, nsub, isub, type, maxsub, ipm,
66
 
     &   ntmp1, ntmp2
 
69
 
 
70
      integer iter, n, n2, nvec, nsub, isub, type, maxsub, 
 
71
     &        ntmp1, ntmp2
67
72
 
68
73
c ... jochen: for convenience, now most arrays have two components.
69
74
c     that might be changed later if memory becomes an issue
70
 
      integer g_y, g_Ay, g_Ax(2), g_r(2), g_r2, g_a, g_bb,
71
 
     &   g_c, g_xold(2), g_Axold(2), g_Ax_im(2)
72
 
      double precision rmax, rmax1, rmax2, acc
 
75
      integer g_y, g_Ay, 
 
76
     &        g_Ax(2), g_r(2), g_r2, 
 
77
     &        g_a, g_bb,g_c, 
 
78
     &        g_xold(2), g_Axold(2), g_Ax_im(2)
 
79
      double precision rmax, rmax1, rmax2, acc,rmx(2)
73
80
      logical converged
74
81
      logical odebug, debug, converge_precond
 
82
 
 
83
      double precision omg(2)
 
84
      integer dsp
75
85
c
76
86
c     =================================================================
77
87
c
78
88
      debug = (.false. .and. ga_nodeid().eq.0) ! for code development
79
89
 
 
90
c      debug=.true.
 
91
 
80
92
c     check input key if we should check for convergence
81
93
c     after the preconditioner has been applied to the residual
82
94
      if (.not. rtdb_get(rtdb, 'aoresponse:precond',    mt_log, 1,
108
120
     &     nvec,CALC_ERR)
109
121
      endif
110
122
      
111
 
c     later we combine the two components to vecors of double
 
123
c     later we combine the two components to vectors of double
112
124
c     length if we have two components, otherwise not:
113
125
      n2 = n
114
126
      if (ncomp.gt.1) n2 = n+n                  
188
200
c     ---------------------
189
201
c     start interation loop
190
202
c     ---------------------
191
 
c     
192
 
      do iter = 1, maxiter
193
 
        
 
203
c   
 
204
c 000000000000 check guess 000000000 START  
 
205
        if (debug) then
 
206
         do ipm=1,ncomp
 
207
             if (ga_nodeid().eq.0)
 
208
     &        write(*,51) iter,ipm
 
209
 51           format('-------g_x-re-guess(',
 
210
     &               i3,',',i3,')--------- START') 
 
211
             call ga_print(g_x(ipm))
 
212
             if (ga_nodeid().eq.0)
 
213
     &        write(*,52) iter,ipm
 
214
 52           format('-------g_x-re-guess(',
 
215
     &               i3,',',i3,')--------- END') 
 
216
             if (lifetime) then
 
217
             if (ga_nodeid().eq.0)
 
218
     &        write(*,53) iter,ipm
 
219
 53           format('-------g_x-im-guess(',
 
220
     &               i3,',',i3,')--------- START') 
 
221
             call ga_print(g_x_im(ipm))
 
222
             if (ga_nodeid().eq.0)
 
223
     &        write(*,54) iter,ipm
 
224
 54           format('-------g_x-im-guess(',
 
225
     &               i3,',',i3,')--------- END') 
 
226
             endif
 
227
         enddo
 
228
        endif ! end-if-debug
 
229
c 000000000000 check guess 000000000 END
 
230
      do iter = 1, maxiter       
194
231
c       
195
232
c ... jochen: here in the iteration loops we keep track
196
233
c       of two components of the solution vector, ipm = 1 and 2
198
235
c       
199
236
        if (odiff) then
200
237
          do ipm = 1,ncomp   
201
 
            call ga_add(1.0d0, g_x(ipm), -1.0d0,
202
 
     &         g_xold(ipm),  g_x(ipm))
 
238
            call ga_add(1.0d0,g_x(ipm), 
 
239
     &                 -1.0d0,g_xold(ipm),
 
240
     &                        g_x(ipm))
203
241
            call ga_sync()
204
 
          enddo
 
242
          enddo ! end-loop-ncomp
205
243
        endif
206
244
c       
207
245
c ... jochen: call product routine with initial or intermediate
208
246
c       solution vector: g_x and g_Ax MUST have two components here
209
247
        
210
248
        if (debug) write (6,*) 'calling product from ga_lkain_2cpl'
211
 
        call product(acc, g_x, g_Ax, g_x_im, g_Ax_im, omega, limag,
212
 
     &     lifetime, gamwidth, ncomp)
 
249
        if (debug) then 
 
250
           if (ga_nodeid().eq.0) then
 
251
             write(*,10) iter
 
252
 10          format('------ g_x-BEF-product(',i3,')---START')
 
253
           endif
 
254
            do ipm=1,ncomp
 
255
             if (ga_nodeid().eq.0)
 
256
     &        write(*,31) iter,ipm
 
257
 31           format('-------g_x-re(',i3,',',i3,')--------- START') 
 
258
             call ga_print(g_x(ipm))
 
259
             if (ga_nodeid().eq.0)
 
260
     &        write(*,32) iter,ipm
 
261
 32           format('-------g_x-re(',i3,',',i3,')--------- END') 
 
262
             if (ga_nodeid().eq.0)
 
263
     &        write(*,33) iter,ipm
 
264
 33           format('-------g_Ax-re(',i3,',',i3,')--------- START') 
 
265
             call ga_print(g_Ax(ipm))
 
266
             if (ga_nodeid().eq.0)
 
267
     &        write(*,34) iter,ipm
 
268
 34           format('-------g_Ax-re(',i3,',',i3,')--------- END') 
 
269
             if (lifetime) then
 
270
             if (ga_nodeid().eq.0)
 
271
     &        write(*,35) iter,ipm
 
272
 35           format('-------g_x-im(',i3,',',i3,')--------- START') 
 
273
              call ga_print(g_x_im(ipm))
 
274
             if (ga_nodeid().eq.0)
 
275
     &        write(*,36) iter,ipm
 
276
 36           format('-------g_x-im(',i3,',',i3,')--------- END') 
 
277
             if (ga_nodeid().eq.0)
 
278
     &        write(*,37) iter,ipm
 
279
 37           format('-------g_Ax-im(',i3,',',i3,')--------- START') 
 
280
              call ga_print(g_Ax_im(ipm))
 
281
             if (ga_nodeid().eq.0)
 
282
     &        write(*,38) iter,ipm
 
283
 38           format('-------g_Ax-im(',i3,',',i3,')--------- END') 
 
284
             endif
 
285
            enddo ! end-loop-ipm
 
286
           if (ga_nodeid().eq.0) then
 
287
             write(*,11) iter
 
288
 11          format('------ g_x-BEF-product(',i3,')---END')
 
289
           endif
 
290
        endif ! end-if-debug
 
291
 
 
292
        call product(acc, 
 
293
     &               g_x   , g_Ax, 
 
294
     &               g_x_im, g_Ax_im, 
 
295
     &               omega, limag,
 
296
     &               lifetime, gamwidth, ncomp)
 
297
        if (debug) then
 
298
           if (ga_nodeid().eq.0) then
 
299
             write(*,12) iter
 
300
 12          format('------ g_x-AFT-product(',i3,')---START')
 
301
           endif
 
302
            do ipm=1,ncomp
 
303
             if (ga_nodeid().eq.0)
 
304
     &        write(*,41) iter,ipm
 
305
 41           format('-------g_x-re(',i3,',',i3,')--------- START') 
 
306
             call ga_print(g_x(ipm))
 
307
             if (ga_nodeid().eq.0)
 
308
     &        write(*,42) iter,ipm
 
309
 42           format('-------g_x-re(',i3,',',i3,')--------- END') 
 
310
             if (ga_nodeid().eq.0)
 
311
     &        write(*,43) iter,ipm
 
312
 43           format('-------g_Ax-re(',i3,',',i3,')--------- START') 
 
313
             call ga_print(g_Ax(ipm))
 
314
             if (ga_nodeid().eq.0)
 
315
     &        write(*,44) iter,ipm
 
316
 44           format('-------g_Ax-re(',i3,',',i3,')--------- END') 
 
317
             if (lifetime) then
 
318
             if (ga_nodeid().eq.0)
 
319
     &        write(*,45) iter,ipm
 
320
 45           format('-------g_x-im(',i3,',',i3,')--------- START') 
 
321
              call ga_print(g_x_im(ipm))
 
322
             if (ga_nodeid().eq.0)
 
323
     &        write(*,46) iter,ipm
 
324
 46           format('-------g_x-im(',i3,',',i3,')--------- END') 
 
325
             if (ga_nodeid().eq.0)
 
326
     &        write(*,47) iter,ipm
 
327
 47           format('-------g_Ax-im(',i3,',',i3,')--------- START') 
 
328
              call ga_print(g_Ax_im(ipm))
 
329
             if (ga_nodeid().eq.0)
 
330
     &        write(*,48) iter,ipm
 
331
 48           format('-------g_Ax-im(',i3,',',i3,')--------- END') 
 
332
             endif
 
333
            enddo ! end-loop-ipm
 
334
           if (ga_nodeid().eq.0) then
 
335
             write(*,13) iter
 
336
 13          format('------ g_x-AFT-product(',i3,')---END')
 
337
           endif
 
338
        endif ! end-if-debug
 
339
 
213
340
        call ga_sync()
214
341
        if (debug) write (6,*) 'returning product from ga_lkain_2cpl'
215
342
 
220
347
        do ipm = 1,ncomp
221
348
          
222
349
          if (odiff) then
223
 
            call ga_add(1.0d0, g_Ax(ipm), 1.0d0,
224
 
     &         g_Axold(ipm), g_Ax(ipm))
225
 
            call ga_add(1.0d0, g_x(ipm),  1.0d0,
226
 
     &         g_xold(ipm),  g_x(ipm))
 
350
            call ga_add(1.0d0, g_Ax(ipm), 
 
351
     &                  1.0d0, g_Axold(ipm), 
 
352
     &                         g_Ax(ipm))
 
353
            call ga_add(1.0d0, g_x(ipm),  
 
354
     &                  1.0d0, g_xold(ipm),  
 
355
     &                         g_x(ipm))
227
356
            call ga_sync()
228
357
            call ga_copy(g_x(ipm), g_xold(ipm))
229
358
            call ga_copy(g_Ax(ipm), g_Axold(ipm))
234
363
c         g_Ax = g_b if the system is solved. During the first cycle,
235
364
c         g_Ax is calculated from the initial guess
236
365
          call ga_add(1.0d0, g_b(ipm),
237
 
     &       -1.0d0, g_Ax(ipm), g_r(ipm)) ! The residual
238
 
          call ga_sync()
239
 
c         
 
366
     &               -1.0d0, g_Ax(ipm), 
 
367
     &                       g_r(ipm)) ! The residual
 
368
          call ga_sync()        
240
369
        enddo                   ! ipm = 1,ncomp
241
370
        call ga_sync()
242
 
 
 
371
          if (debug) then
 
372
           do ipm=1,ncomp
 
373
            if (ga_nodeid().eq.0) then
 
374
             write(*,49) ipm,iter
 
375
 49          format('----g_Ax(',i3,',',i3,')-----START')
 
376
            endif
 
377
            call ga_print(g_Ax(ipm))
 
378
            if (ga_nodeid().eq.0) then
 
379
             write(*,50) ipm,iter
 
380
 50          format('----g_Ax(',i3,',',i3,')-----END')
 
381
            endif
 
382
            if (ga_nodeid().eq.0) then
 
383
             write(*,2779) ipm,iter
 
384
 2779        format('----g_b(',i3,',',i3,')-----START')
 
385
            endif
 
386
            call ga_print(g_b(ipm))
 
387
            if (ga_nodeid().eq.0) then
 
388
             write(*,2880) ipm,iter
 
389
 2880        format('----g_b(',i3,',',i3,')-----END')
 
390
            endif
 
391
            if (ga_nodeid().eq.0) then
 
392
             write(*,2782) iter,ipm
 
393
 2782        format('----g_r(',i3,',',i3,')-----START')
 
394
            endif
 
395
            call ga_print(g_r(ipm))
 
396
            if (ga_nodeid().eq.0) then
 
397
             write(*,2783) iter,ipm
 
398
 2783        format('----g_r(',i3,',',i3,')-----END')
 
399
            endif
 
400
           enddo ! end-loop-ipm
 
401
          endif ! end-if-debug      
243
402
c       convergence checking:
244
403
c       find the largest element of the residual either 
245
404
c       before or after the call to the preconditioner
246
 
 
247
 
        if (converge_precond) then
248
 
          call precond(g_r(1),  -omega)
249
 
          if (ncomp.gt.1) then
250
 
            call precond(g_r(2),   omega)
251
 
          endif
252
 
        endif
253
 
 
254
 
        call ga_maxelt(g_r(1), rmax1)
255
 
        if (ncomp.gt.1) then
256
 
          call ga_maxelt(g_r(2), rmax2)
257
 
        else
258
 
          rmax2 = 0d0
259
 
        endif
260
 
        rmax = max(rmax1, rmax2)  
 
405
         omg(1)=-omega
 
406
         omg(2)= omega
 
407
 
 
408
         if (converge_precond) then
 
409
          do ipm=1,ncomp
 
410
           call precond(g_r(ipm),omg(ipm))
 
411
          enddo ! end-loop-ipm
 
412
         endif ! end-if-converge_precond
 
413
         rmx(1)=0.0d0
 
414
         rmx(2)=0.0d0
 
415
          do ipm=1,ncomp
 
416
           call ga_maxelt(g_r(ipm),rmx(ipm))   
 
417
          enddo ! end-loop-ipm
 
418
 
 
419
         rmax1=rmx(1)
 
420
         rmax2=rmx(2)
 
421
 
 
422
         rmax = max(rmax1, rmax2)  
261
423
      
262
424
        if (oprint .and. ga_nodeid().eq.0) then
263
425
          write(6,3) iter, nsub+nvec, rmax, util_wallsec()
275
437
c       (there were only two in ga_lkain, one with g_aX, one with g_r)
276
438
c       for array g_r the preconditioner call is only necessary
277
439
c       in case converge_precond is .false.
278
 
 
279
 
        call precond(g_Ax(1), -omega)
280
 
        if (.not.converge_precond) call precond(g_r(1), -omega)
281
 
        if (ncomp.gt.1) then
282
 
          call precond(g_Ax(2),  omega)
283
 
          if (.not.converge_precond) call precond(g_r(2), omega)
284
 
        endif       
 
440
        omg(1)=-omega
 
441
        omg(2)= omega
 
442
        do ipm=1,ncomp
 
443
         call precond(g_Ax(ipm),omg(ipm))
 
444
         if (.not.converge_precond) 
 
445
     &   call precond(g_r(ipm) ,omg(ipm))
 
446
        enddo ! end-loop-ipm
 
447
  
285
448
        call ga_sync()
286
449
        
287
 
c       Copy the vectors to the subspace work area
288
 
c       
289
 
        call ga_copy_patch('n', 
290
 
     $     g_Ax(1), 1, n, 1, nvec, 
291
 
     $     g_Ay, 1, n, nsub+1, nsub+nvec)
292
 
        call ga_sync()
293
 
        if (ncomp.gt.1) then
294
 
          call ga_copy_patch('n', 
295
 
     $       g_Ax(2), 1, n, 1, nvec, 
296
 
     $       g_Ay, n+1, n2, nsub+1, nsub+nvec)
297
 
        endif
298
 
        call ga_copy_patch('n', 
299
 
     $     g_x(1), 1, n, 1, nvec, 
300
 
     $     g_y, 1, n, nsub+1, nsub+nvec)
301
 
        call ga_sync()
302
 
        if (ncomp.gt.1) then
303
 
          call ga_copy_patch('n', 
304
 
     $       g_x(2), 1, n, 1, nvec, 
305
 
     $       g_y, n+1, n2, nsub+1, nsub+nvec)
306
 
        endif
307
 
 
308
 
c       g_r2 is needed below for multiplication
309
 
        call ga_copy_patch('n', 
310
 
     $     g_r(1), 1, n, 1, nvec, 
311
 
     $     g_r2, 1, n, 1, nvec)
312
 
        if (ncomp.gt.1) then
313
 
          call ga_copy_patch('n', 
314
 
     $       g_r(2), 1, n, 1, nvec, 
315
 
     $       g_r2, n+1, n2, 1, nvec)
316
 
        endif
 
450
c       Copy the vectors to the subspace work area      
 
451
        do ipm=1,ncomp
 
452
         dsp=n*(ipm-1)
 
453
         call ga_copy_patch('n', 
 
454
     $                      g_Ax(ipm),1    ,n    ,1     ,nvec, 
 
455
     $                      g_Ay     ,1+dsp,n+dsp,nsub+1,nsub+nvec)
 
456
         call ga_copy_patch('n', 
 
457
     $                      g_x(ipm),1    ,n    ,1     ,nvec, 
 
458
     $                      g_y     ,1+dsp,n+dsp,nsub+1,nsub+nvec)
 
459
         call ga_copy_patch('n', 
 
460
     $                      g_r(ipm),1    ,n    ,1,nvec, 
 
461
     $                      g_r2    ,1+dsp,n+dsp,1,nvec)
 
462
        enddo ! end-loop-ipm
317
463
        
318
464
        nsub = nsub + nvec
319
465
c       
334
480
        call ga_zero(g_c)
335
481
        call ga_sync()
336
482
        call ga_dgemm('t','n',nsub,nsub,n2,1.0d0,
337
 
     &     g_y,g_Ay,0.0d0,g_a)
 
483
     &                g_y,g_Ay,0.0d0,g_a)
338
484
        call ga_dgemm('t','n',nsub,nvec,n2,1.0d0,
339
 
     &     g_y,g_r2,0.0d0,g_bb)
 
485
     &                g_y,g_r2,0.0d0,g_bb)
340
486
        call ga_sync()
341
487
        if (odebug) call ga_print(g_a)
342
488
        if (odebug) call ga_print(g_c)
357
503
c       
358
504
        call ga_sync()
359
505
        call ga_dgemm('n','n',n2,nvec,nsub,-1.0d0,
360
 
     &     g_Ay,g_c,1.0d0,g_r2)
 
506
     &                g_Ay,g_c,1.0d0,g_r2)
361
507
        if (odebug) then
362
508
          write(6,*) ' The update in the complement '
363
509
          call ga_print(g_r2)
364
510
        end if
365
511
c       
366
512
c       copy components of g_r2 into g_r before adding g_r to  g_x
367
 
c
368
 
        call ga_sync()
369
 
        call ga_copy_patch('n', 
370
 
     $     g_r2, 1, n, 1, nvec, 
371
 
     $     g_r(1), 1, n, 1,nvec)
372
 
        if (ncomp.gt.1) then
373
 
          call ga_copy_patch('n', 
374
 
     $       g_r2, n+1, n2, 1, nvec, 
375
 
     $       g_r(2), 1, n, 1,nvec)
376
 
        endif
377
 
        call ga_sync()
378
 
c       
 
513
        do ipm=1,ncomp
 
514
         dsp=n*(ipm-1)
 
515
         call ga_copy_patch('n', 
 
516
     $                     g_r2    ,1+dsp,n+dsp,1,nvec, 
 
517
     $                     g_r(ipm),1    ,n    ,1,nvec)      
 
518
        enddo ! end-loop-ipm
 
519
        call ga_sync()       
379
520
        do ipm = 1,ncomp
380
 
          call ga_add(1.0d0, g_r(ipm), 1.0d0, g_x(ipm), g_x(ipm))
381
 
        enddo
 
521
          call ga_add(1.0d0, g_r(ipm), 
 
522
     &                1.0d0, g_x(ipm), 
 
523
     &                       g_x(ipm))
 
524
        enddo ! end-loop-ipm
382
525
        
383
526
        call ga_sync()
384
527
        call ga_dgemm('n','n',n2,nvec,nsub,1.0d0,
385
 
     &     g_y,g_c,0.0d0,g_r2)
 
528
     &                g_y,g_c,0.0d0,g_r2)
 
529
 
386
530
        call ga_sync()
387
531
        if (odebug) then
388
532
          write(6,*) ' The update in the subspace '
389
533
          call ga_print(g_r2)
390
534
        end if
391
535
c       
392
 
c       copy components of g_r2 into g_r before adding g_r to  g_x
393
 
c
394
 
        call ga_copy_patch('n', 
395
 
     $     g_r2, 1, n, 1, nvec, 
396
 
     $     g_r(1), 1, n, 1,nvec)
397
 
        if (ncomp.gt.1) then
398
 
          call ga_copy_patch('n', 
399
 
     $       g_r2, n+1, n2, 1, nvec, 
400
 
     $       g_r(2), 1, n, 1,nvec)
401
 
        endif
 
536
c       copy components of g_r2 into g_r before adding g_r to  g_x    
 
537
        do ipm=1,ncomp
 
538
         dsp=n*(ipm-1)
 
539
         call ga_copy_patch('n', 
 
540
     $                     g_r2    ,1+dsp,n+dsp,1,nvec, 
 
541
     $                     g_r(ipm),1    ,n    ,1,nvec)      
 
542
        enddo ! end-loop-ipm
402
543
        call ga_sync()
403
544
        do ipm = 1,ncomp
404
 
          call ga_add(1.0d0, g_r(ipm), 1.0d0, g_x(ipm), g_x(ipm))
405
 
        enddo
 
545
          call ga_add(1.0d0, g_r(ipm), 
 
546
     &                1.0d0, g_x(ipm), 
 
547
     &                       g_x(ipm))
 
548
        enddo ! end-loop-ipm
406
549
        call ga_sync()
407
550
c        
408
551
        if (.not. ga_destroy(g_a)) call errquit
507
650
c     
508
651
      end
509
652
 
510
 
c     ******************************************************************
511
 
c     ******************************************************************
512
 
 
513
 
      subroutine ga_lkain_2cpl3_damp(rtdb, g_x, g_b, g_x_im, g_b_im,
514
 
     &   product, precond, 
515
 
     $   tol, mmaxsub, maxiter, odiff, oprint, omega, limag,
516
 
     &   lifetime, gamwidth, ncomp)
 
653
      subroutine ga_lkain_2cpl3_damp(rtdb, 
 
654
     &                               g_x, g_b, 
 
655
     &                               g_x_im, g_b_im,
 
656
     &                               product, precond, 
 
657
     $                               tol, mmaxsub, maxiter, 
 
658
     &                               odiff, oprint, omega, limag,
 
659
     &                               lifetime, gamwidth, ncomp)
 
660
c  Written by J. Autschbach, SUNY Buffalo
 
661
c  Improvements made
 
662
c          by F. Aquino,     Northwestern University 
 
663
c          03-15-12
 
664
c  Note.- Modifying/Improving ga_lkain_2cpl3_damp() as of old src code
 
665
c         appeared on download src 11-20-12
 
666
c --> Experimental (not published yet)
517
667
      implicit none
518
668
#include "errquit.fh"
519
669
#include "mafdecls.fh"
587
737
      double precision rmax, rmax1, rmax2, acc
588
738
      logical converged
589
739
      logical odebug, debug, converge_precond
 
740
      double complex val_cmplx
 
741
      logical debug1
 
742
      integer p1,p2,m1,m2,stat_solve
 
743
      
 
744
c      integer l_zre,k_zre,g_z,
 
745
c     &        l_zim,k_zim
 
746
      external copy_rtor2,copy_r2tor,copy_AxxtoAyy,
 
747
     &         update_g_x1
 
748
 
 
749
      debug1 = .false. ! FA no  printouts
 
750
c      debug1 = .true.  ! FA yes printouts
590
751
c
591
752
c     =================================================================
592
 
c
 
753
 
593
754
      debug = (.false. .and. ga_nodeid().eq.0) ! for code development
594
755
 
595
756
c     check input key if we should check for convergence
623
784
 
624
785
      call ga_inquire(g_x(1), type, n, nvec)
625
786
 
 
787
c ------- create (zre,zim) ---------- START
 
788
c      if (.not.MA_Push_Get(mt_dbl,n,'hessv jfacs',l_zre,k_zre))
 
789
c     &     call errquit('sh-fockbld2: cannot allocate jfac',
 
790
c     &                  n, MA_ERR)
 
791
c      if (.not.MA_Push_Get(mt_dbl,n,'hessv kfacs',l_zim,k_zim))
 
792
c     &     call errquit('sh-fockbld2: cannot allocate kfac',
 
793
c     &                  n, MA_ERR)
 
794
c        if (.not. ga_create(MT_DCPL,n,nvec, 'lkain_2cpl: z',
 
795
c     $     0, 0, g_z))
 
796
c     $     call errquit('lkain: failed allocating z', nvec,
 
797
c     &     GA_ERR)
 
798
c           call ga_zero(g_z)
 
799
c ------- create (zre,zim) ---------- END
626
800
      if (ncomp.gt.1) then
627
801
        call ga_inquire(g_x(2), type, ntmp1, ntmp2)
628
802
        
636
810
c     length and combine again Re and Im, i.e. 
637
811
c     the dimension is up to 4*n
638
812
 
 
813
c ========= to be removed ============ start
639
814
      n2 = n
640
815
      if (ncomp.gt.1) n2 = 2 * n   
641
816
      n3 = n
645
820
      if (ncomp.gt.1 .or. lifetime) n4 = 2* n
646
821
      if (lifetime .and. ncomp.gt.1) n4 = 4 * n
647
822
      if (debug) write (6,*) 'n1n2n3n4',n,n2,n3,n4
 
823
c ========= to be removed ============ end
648
824
 
649
825
      maxsub = mmaxsub          ! So don't modify input scalar arg
650
826
      if (maxsub .lt. 3*nvec) maxsub = 3*nvec
749
925
c       
750
926
        if (odiff) then
751
927
          do ipm = 1,ncomp   
752
 
            call ga_add(1.0d0, g_x(ipm), -1.0d0,
753
 
     &         g_xold(ipm),  g_x(ipm))
 
928
            call ga_add( 1.0d0,g_x(ipm), 
 
929
     &                  -1.0d0,g_xold(ipm),  
 
930
     &                         g_x(ipm))
754
931
            call ga_sync()
755
932
          enddo
756
933
        endif
758
935
c ... jochen: call product routine with initial or intermediate
759
936
c       solution vector: g_x and g_Ax MUST have dimension two here
760
937
c       even if only one of them is used
 
938
          if (debug) then   
 
939
            do ipm=1,ncomp
 
940
             if (ga_nodeid().eq.0) then
 
941
              write(*,112) iter,ipm
 
942
 112          format('------ prod-g_x-1(',i3,',',i3,')------ START')
 
943
             endif
 
944
             call ga_print(g_x(ipm))
 
945
             if (ga_nodeid().eq.0) then
 
946
              write(*,113) iter,ipm
 
947
 113          format('------ prod-g_x-1(',i3,',',i3,')------ END')
 
948
             endif
 
949
             if (ga_nodeid().eq.0) then
 
950
              write(*,114) iter,ipm
 
951
 114          format('------ prod-g_Ax-1(',i3,',',i3,')------ START')
 
952
             endif
 
953
             call ga_print(g_Ax(ipm))
 
954
             if (ga_nodeid().eq.0) then
 
955
              write(*,115) iter,ipm
 
956
 115          format('------ prod-g_Ax-1(',i3,',',i3,')------ END')
 
957
             endif
 
958
             if (lifetime) then
 
959
              call ga_print(g_x_im(ipm))
 
960
              call ga_print(g_Ax_im(ipm))
 
961
             endif
 
962
            enddo ! end-loop-ipm
 
963
         endif ! end-if-debug
761
964
        
762
965
        if (debug) write (6,*)
763
966
     &     'calling product from ga_lkain_2cpl_damp'
764
967
        call product(acc, g_x, g_Ax,  g_x_im, g_Ax_im, omega, limag,
765
968
     &     lifetime, gamwidth, ncomp)
 
969
         if (debug) then
 
970
            do ipm=1,ncomp
 
971
             if (ga_nodeid().eq.0) then
 
972
              write(*,116) iter,ipm
 
973
 116          format('------ prod-g_x-2(',i3,',',i3,')------ START')
 
974
             endif
 
975
             call ga_print(g_x(ipm))
 
976
             if (ga_nodeid().eq.0) then
 
977
              write(*,117) iter,ipm
 
978
 117          format('------ prod-g_x-2(',i3,',',i3,')------ END')
 
979
             endif
 
980
             if (ga_nodeid().eq.0) then
 
981
              write(*,118) iter,ipm
 
982
 118          format('------ prod-g_Ax-2(',i3,',',i3,')------ START')
 
983
             endif
 
984
             call ga_print(g_Ax(ipm))
 
985
             if (ga_nodeid().eq.0) then
 
986
              write(*,119) iter,ipm
 
987
 119          format('------ prod-g_Ax-2(',i3,',',i3,')------ END')
 
988
             endif
 
989
             enddo ! end-loop-ipm
 
990
           endif ! end-if-debug
 
991
c             if (iter.eq.3) then
 
992
c               if (ga_nodeid().eq.0)
 
993
c     &          write(*,*) 'FA-STOP-convergence at 3rd iter'
 
994
c               stop
 
995
c             endif              
 
996
 
766
997
        if (debug) write (6,*)
767
998
     &     'returning product from ga_lkain_2cpl_damp'
768
999
 
775
1006
          
776
1007
          if (odiff) then
777
1008
c           jochen: odiff stuff presently ignored
778
 
            call ga_add(1.0d0, g_Ax(ipm), 1.0d0,
779
 
     &         g_Axold(ipm), g_Ax(ipm))
780
 
            call ga_add(1.0d0, g_x(ipm),  1.0d0,
781
 
     &         g_xold(ipm),  g_x(ipm))
 
1009
            call ga_add(1.0d0,g_Ax(ipm), 
 
1010
     &                  1.0d0,g_Axold(ipm), 
 
1011
     &                        g_Ax(ipm))
 
1012
            call ga_add(1.0d0,g_x(ipm),  
 
1013
     &                  1.0d0,g_xold(ipm),
 
1014
     &                        g_x(ipm))
782
1015
            call ga_sync()
783
1016
            call ga_copy(g_x(ipm), g_xold(ipm))
784
1017
            call ga_copy(g_Ax(ipm), g_Axold(ipm))
798
1031
c         and assigned real and imaginary parts accordingly)
799
1032
 
800
1033
          call ga_sync()
801
 
          call ga_add(1.0d0, g_b(ipm),
802
 
     &       -1.0d0, g_Ax(ipm), g_r(ipm)) ! The residual, Real part
 
1034
c FA: Step 1:
 
1035
          call ga_add( 1.0d0,g_b(ipm),
 
1036
     &                -1.0d0,g_Ax(ipm), 
 
1037
     &                       g_r(ipm)) ! The residual, Real part
 
1038
 
 
1039
          if (debug) then
 
1040
             if (ga_nodeid().eq.0) then
 
1041
              write(*,120) iter,ipm
 
1042
 120          format('------ prod-g_b(',i3,',',i3,')------ START')
 
1043
             endif
 
1044
             call ga_print(g_b(ipm))
 
1045
             if (ga_nodeid().eq.0) then
 
1046
              write(*,121) iter,ipm
 
1047
 121          format('------ prod-g_b(',i3,',',i3,')------ END')
 
1048
             endif
 
1049
             if (ga_nodeid().eq.0) then
 
1050
              write(*,122) iter,ipm
 
1051
 122          format('------ prod-g_r(',i3,',',i3,')------ START')
 
1052
             endif
 
1053
             call ga_print(g_r(ipm))
 
1054
             if (ga_nodeid().eq.0) then
 
1055
              write(*,123) iter,ipm
 
1056
 123          format('------ prod-g_r(',i3,',',i3,')------ END')
 
1057
             endif
 
1058
          endif ! end-if-debug
803
1059
 
804
1060
          if (lifetime) then
805
 
            call ga_add(1.0d0, g_b_im(ipm),
806
 
     &         -1.0d0, g_Ax_im(ipm), g_r_im(ipm)) ! The residual, Im part
 
1061
            call ga_add( 1.0d0,g_b_im(ipm),
 
1062
     &                  -1.0d0,g_Ax_im(ipm), 
 
1063
     &                         g_r_im(ipm)) ! The residual, Im part
807
1064
          endif
808
1065
c         
809
1066
        enddo                   ! ipm = 1,ncomp
810
 
        call ga_sync()
811
1067
 
812
1068
c       convergence checking:
813
1069
c       find the largest element of the residual either 
814
1070
c       before or after the call to the preconditioner
815
1071
 
816
1072
        if (converge_precond) then
817
 
          call precond(g_r(1),  g_r_im(1),  -omega, gamwidth)
818
 
          if (ncomp.gt.1) then
819
 
            call precond(g_r(2),  g_r_im(2),   omega, gamwidth)
820
 
          endif          
 
1073
         do ipm=1,ncomp
 
1074
          call precond(g_r(ipm),g_r_im(ipm),-omega,gamwidth)
 
1075
         enddo ! end-loop-ipm
821
1076
        endif
822
 
        call ga_sync()
823
 
c
 
1077
 
824
1078
        call ga_maxelt(g_r(1), rmax1)
825
1079
        if (ncomp.gt.1) then
826
1080
          call ga_maxelt(g_r(2), rmax2)
863
1117
          if (.not.converge_precond) 
864
1118
     &       call precond(g_r(2),  g_r_im(2),   omega, gamwidth)
865
1119
        endif
866
 
        call ga_sync()
 
1120
 
867
1121
        if (debug) write (6,*) 'lkain3_damp: back from precond'
868
1122
                
869
1123
c       Copy the vectors to the subspace work area
870
 
c       
871
 
        call ga_copy_patch('n', 
872
 
     $     g_Ax(1), 1, n, 1, nvec, 
873
 
     $     g_Ay, 1, n, nsub+1, nsub+nvec)
874
 
        if (ncomp.gt.1) call ga_copy_patch('n', 
875
 
     $     g_Ax(2), 1, n, 1, nvec, 
876
 
     $     g_Ay, n+1, n2, nsub+1, nsub+nvec)
877
 
        if (lifetime) then
878
 
          call ga_copy_patch('n', 
879
 
     $       g_Ax_im(1), 1, n, 1, nvec, 
880
 
     $       g_Ay, n2+1, n3, nsub+1, nsub+nvec)
881
 
          if (ncomp.gt.1) call ga_copy_patch('n', 
882
 
     $       g_Ax_im(2), 1, n, 1, nvec, 
883
 
     $       g_Ay, n3+1, n4, nsub+1, nsub+nvec)
884
 
        endif                   ! lifetime
885
 
        call ga_copy_patch('n', 
886
 
     $     g_x(1), 1, n, 1, nvec, 
887
 
     $     g_y, 1, n, nsub+1, nsub+nvec)
888
 
        if (ncomp.gt.1)  call ga_copy_patch('n', 
889
 
     $     g_x(2), 1, n, 1, nvec, 
890
 
     $     g_y, n+1, n2, nsub+1, nsub+nvec)
891
 
        if (lifetime) then
892
 
          call ga_copy_patch('n', 
893
 
     $       g_x_im(1), 1, n, 1, nvec, 
894
 
     $       g_y, n2+1, n3, nsub+1, nsub+nvec)
895
 
          if (ncomp.gt.1) call ga_copy_patch('n', 
896
 
     $       g_x_im(2), 1, n, 1, nvec, 
897
 
     $       g_y, n3+1, n4, nsub+1, nsub+nvec)
898
 
        endif                   ! lifetime
899
 
        call ga_sync()
900
 
        if (debug) write (6,*) 'lkain3_damp: vec to sub complete'
901
 
 
902
 
c       g_r2 is needed below for multiplication
903
 
        call ga_copy_patch('n', 
904
 
     $     g_r(1), 1, n, 1, nvec, 
905
 
     $     g_r2, 1, n, 1, nvec)
906
 
        if (ncomp.gt.1) call ga_copy_patch('n', 
907
 
     $     g_r(2), 1, n, 1, nvec, 
908
 
     $     g_r2, n+1, n2, 1, nvec)
909
 
        if (lifetime) then
910
 
          call ga_copy_patch('n', 
911
 
     $       g_r_im(1), 1, n, 1, nvec, 
912
 
     $       g_r2, n2+1, n3, 1, nvec)
913
 
          if (ncomp.gt.1) call ga_copy_patch('n', 
914
 
     $       g_r_im(2), 1, n, 1, nvec, 
915
 
     $       g_r2, n3+1, n4, 1, nvec)
916
 
        endif                   ! lifetime
917
 
        call ga_sync()
918
 
        if (debug) write (6,*) 'lkain3_damp: r2 complete'
919
 
        
 
1124
c   
 
1125
c ---- FA-copy ((Ax,Ax_im),(x,x_im)) -> (Ay,y) ------ START
 
1126
c FA: Step 2:
 
1127
        if (debug1) then
 
1128
        if (ga_nodeid().eq.0)
 
1129
     &   write(*,*) '---------g_x-0(',iter,')-----START' 
 
1130
        call ga_print(g_x(1))
 
1131
        if (ga_nodeid().eq.0)
 
1132
     &   write(*,*) '---------g_x-0(',iter,')-----END'    
 
1133
        if (ga_nodeid().eq.0)
 
1134
     &   write(*,*) '---------g_x_im-0(',iter,')-----START' 
 
1135
        call ga_print(g_x_im(1))
 
1136
        if (ga_nodeid().eq.0)
 
1137
     &   write(*,*) '---------g_x_im-0(',iter,')-----END'        
 
1138
        if (ga_nodeid().eq.0)
 
1139
     &   write(*,*) '---------g_Ax-0(',iter,')-----START' 
 
1140
        call ga_print(g_Ax(1))
 
1141
        if (ga_nodeid().eq.0)
 
1142
     &   write(*,*) '---------g_Ax-0(',iter,')-----END'          
 
1143
        if (ga_nodeid().eq.0)
 
1144
     &   write(*,*) '---------g_r-0(',iter,')-----START' 
 
1145
        call ga_print(g_r(1))
 
1146
        if (ga_nodeid().eq.0)
 
1147
     &   write(*,*) '---------g_r-0(',iter,')-----END'  
 
1148
        endif ! end-if-debug1  
 
1149
     
 
1150
        call copy_AxxtoAyy(g_Ax,g_Ax_im,
 
1151
     &                     g_x,g_x_im,
 
1152
     &                     g_Ay,g_y,
 
1153
     &                     nvec,
 
1154
     &                     ncomp,
 
1155
     &                     nsub,
 
1156
     &                     n,
 
1157
     &                     lifetime)
 
1158
 
 
1159
c ---- FA-copy ((Ax,Ax_im),(x,x_im)) -> (Ay,y) ------ END
 
1160
c ----- FA--- (g_r,g_r_im) -> g_r2 ---------- START
 
1161
c FA: Step 3:
 
1162
         call copy_rtor2(g_r2,
 
1163
     &                   g_r,
 
1164
     &                   g_r_im,
 
1165
     &                   ncomp,
 
1166
     &                   nvec,
 
1167
     &                   n,
 
1168
     &                   lifetime)
 
1169
c ----- FA--- (g_r,g_r_im) -> g_r2 ---------- END 
 
1170
        if (debug1) then
 
1171
        if (ga_nodeid().eq.0)
 
1172
     &   write(*,*) '---------g_y-0(',iter,')-----START' 
 
1173
        call ga_print(g_y)
 
1174
        if (ga_nodeid().eq.0)
 
1175
     &   write(*,*) '---------g_y-0(',iter,')-----END'          
 
1176
        if (ga_nodeid().eq.0)
 
1177
     &   write(*,*) '---------g_Ay-0(',iter,')-----START' 
 
1178
        call ga_print(g_Ay)
 
1179
        if (ga_nodeid().eq.0)
 
1180
     &   write(*,*) '---------g_Ay-0(',iter,')-----END'          
 
1181
        if (ga_nodeid().eq.0)
 
1182
     &   write(*,*) '---------g_r2-0(',iter,')-----START' 
 
1183
        call ga_print(g_r2)
 
1184
        if (ga_nodeid().eq.0)
 
1185
     &   write(*,*) '---------g_r2-0(',iter,')-----END'    
 
1186
        endif ! end-if-debug1
 
1187
 
920
1188
        nsub = nsub + nvec
921
1189
c       
922
1190
c       Form and solve the subspace equations using SVD in order
934
1202
        call ga_zero(g_a)
935
1203
        call ga_zero(g_bb)
936
1204
        call ga_zero(g_c)
937
 
        call ga_sync()
938
1205
        call ga_dgemm('t','n',nsub,nsub,n4,1.0d0,
939
 
     &     g_y,g_Ay,0.0d0,g_a)
940
 
        call ga_sync()
 
1206
     &                g_y,g_Ay,0.0d0,g_a)
941
1207
        call ga_dgemm('t','n',nsub,nvec,n4,1.0d0,
942
 
     &     g_y,g_r2,0.0d0,g_bb)
943
 
        call ga_sync()
 
1208
     &                g_y,g_r2,0.0d0,g_bb)
 
1209
 
944
1210
        if (odebug) call ga_print(g_a)
945
1211
        if (odebug) call ga_print(g_c)
946
1212
c       
952
1218
c       realistic examples (maxsub << n) there is only a little
953
1219
c       advantage and in the precence of real noise in the products
954
1220
c       screening with a realistic threshold is important.
955
 
c       
 
1221
        if (debug1) then
 
1222
        if (ga_nodeid().eq.0)
 
1223
     &   write(*,*) '---------g_a(',iter,')-----START' 
 
1224
        call ga_print(g_a)
 
1225
        if (ga_nodeid().eq.0)
 
1226
     &   write(*,*) '---------g_a(',iter,')-----END'
 
1227
        if (ga_nodeid().eq.0)
 
1228
     &   write(*,*) '---------g_bb(',iter,')-----START' 
 
1229
        call ga_print(g_bb)
 
1230
        if (ga_nodeid().eq.0)
 
1231
     &   write(*,*) '---------g_bb(',iter,')-----END'
 
1232
        endif ! end-if-debug1
 
1233
c ++++++++++++++++++++++++++++++++++++++++++++++++++++    
 
1234
c ++++++++++++++ real linear solver ++++++++++++++++++ 
956
1235
        call ga_svd_solve_seq(g_a,g_bb,g_c,1d-14)
 
1236
 
 
1237
c        if (ga_nodeid().eq.0)
 
1238
c     &   write(*,*) 'FA-test ga_lu_solve ...'
 
1239
c        call ga_lu_solve('n',g_a,g_bb)
 
1240
 
 
1241
c        stat_solve=ga_solve(g_a,g_bb)
 
1242
c        if (ga_nodeid().eq.0) then
 
1243
c         write(*,14) stat_solve
 
1244
c 14      format('FA-test ga_lu_solve(Choleski->LU) stat=',i3)
 
1245
c        endif
 
1246
c        call ga_copy(g_bb,g_c)
 
1247
c ++++++++++++++++++++++++++++++++++++++++++++++++++++
 
1248
c ++++++++++++++++++++++++++++++++++++++++++++++++++++
 
1249
        if (debug1) then
 
1250
        if (ga_nodeid().eq.0)
 
1251
     &   write(*,*) '---------g_c(',iter,')-----START' 
 
1252
        call ga_print(g_c)
 
1253
        if (ga_nodeid().eq.0)
 
1254
     &   write(*,*) '---------g_c(',iter,')-----END'
 
1255
        endif ! end-if-debug1
 
1256
        
957
1257
        if (odebug) call ga_print(g_c)
958
1258
c       
959
 
c       Form and add the correction, in parts, onto the solution
960
 
c       
961
 
        call ga_sync()
 
1259
c       Form and add the correction, in parts, onto the solution      
 
1260
c FA: Step 5:
 
1261
        if (debug1) then
 
1262
        if (ga_nodeid().eq.0)
 
1263
     &   write(*,*) '---------g_r2-BEF(',iter,')-----START' 
 
1264
        call ga_print(g_r2)
 
1265
        if (ga_nodeid().eq.0)
 
1266
     &   write(*,*) '---------g_r2-BEF(',iter,')-----END'
 
1267
        endif ! end-if-debug1
 
1268
 
962
1269
        call ga_dgemm('n','n',n4,nvec,nsub,-1.0d0,
963
 
     &     g_Ay,g_c,1.0d0,g_r2)
964
 
        call ga_sync()
965
 
        if (odebug) then
966
 
          write(6,*) ' The update in the complement '
967
 
          call ga_print(g_r2)
968
 
        end if
969
 
        if (debug) write (6,*) 'lkain3_damp: after dgemm'
 
1270
     &                g_Ay,g_c,1.0d0,g_r2)
 
1271
 
 
1272
        if (debug1) then
 
1273
        if (ga_nodeid().eq.0)
 
1274
     &   write(*,*) '---------g_r2-AFT(',iter,')-----START' 
 
1275
        call ga_print(g_r2)
 
1276
        if (ga_nodeid().eq.0)
 
1277
     &   write(*,*) '---------g_r2-AFT(',iter,')-----END'
 
1278
        if (ga_nodeid().eq.0)
 
1279
     &   write(*,*) '---------g_x-BEF-1(',iter,')-----START' 
 
1280
        call ga_print(g_x)
 
1281
        if (ga_nodeid().eq.0)
 
1282
     &   write(*,*) '---------g_x-BEF-1(',iter,')-----END'
 
1283
        endif ! end-if-debug1
970
1284
c       
971
1285
c       copy components of g_r2 into g_r before adding g_r to  g_x
972
 
        call ga_copy_patch('n', 
973
 
     $     g_r2, 1, n, 1, nvec, 
974
 
     $     g_r(1), 1, n, 1,nvec)
975
 
        if (ncomp.gt.1) call ga_copy_patch('n', 
976
 
     $     g_r2, n+1, n2, 1, nvec, 
977
 
     $     g_r(2), 1, n, 1,nvec)
978
 
        if (lifetime) then
979
 
          call ga_copy_patch('n', 
980
 
     $       g_r2, n2+1, n3, 1, nvec, 
981
 
     $       g_r_im(1), 1, n, 1,nvec)
982
 
          if (ncomp.gt.1) call ga_copy_patch('n', 
983
 
     $       g_r2, n3+1, n4, 1, nvec, 
984
 
     $       g_r_im(2), 1, n, 1,nvec)
985
 
        endif                   ! lifetime
986
 
        call ga_sync()
987
 
        if (debug) write (6,*) 'lkain3_damp: r2 copied'
988
 
c               
989
 
        do ipm = 1,ncomp
990
 
          call ga_add(1.0d0, g_r(ipm), 1.0d0, g_x(ipm), g_x(ipm))
991
 
          if (lifetime)
992
 
     &       call ga_add(1.0d0, g_r_im(ipm),
993
 
     &       1.0d0, g_x_im(ipm), g_x_im(ipm))
994
 
        enddo
995
 
        call ga_sync()
996
 
c       
 
1286
       call update_g_x1(g_r2,    ! in   :
 
1287
     &                 g_x,     ! in/ou: updated
 
1288
     &                 g_x_im,  ! in/ou: updated
 
1289
     &                 ncomp,   ! in   : nr. components
 
1290
     &                 nvec,    ! in   : (x,y,z) = 3
 
1291
     &                 n,       ! in   : size of block (e.g. RE1,RE2,IM1 or IM2)
 
1292
     &                 lifetime)! in   : = .true. if complex
 
1293
        if (debug1) then
 
1294
        if (ga_nodeid().eq.0)
 
1295
     &   write(*,*) '---------g_x-AFT-1(',iter,')-----START' 
 
1296
        call ga_print(g_x)
 
1297
        if (ga_nodeid().eq.0)
 
1298
     &   write(*,*) '---------g_x-AFT-1(',iter,')-----END'
 
1299
        endif ! end-if-debug1
 
1300
c FA: Step 8:
997
1301
        call ga_dgemm('n','n',n4,nvec,nsub,1.0d0,
998
 
     &     g_y,g_c,0.0d0,g_r2)
999
 
        call ga_sync()
1000
 
        if (odebug) then
1001
 
          write(6,*) ' The update in the subspace '
1002
 
          call ga_print(g_r2)
1003
 
        end if
1004
 
        if (debug) write (6,*) 'lkain3_damp: subspace updated'
1005
 
        
1006
 
c       
 
1302
     &                g_y,g_c,0.0d0,g_r2)
 
1303
        if (debug1) then
 
1304
        if (ga_nodeid().eq.0)
 
1305
     &   write(*,*) '---------g_y c(',iter,')-----START' 
 
1306
        call ga_print(g_r2)
 
1307
        if (ga_nodeid().eq.0)
 
1308
     &   write(*,*) '---------g_y c(',iter,')-----END'
 
1309
        endif ! end-if-debug1
1007
1310
c       copy components of g_r2 into g_r before adding g_r to  g_x
1008
 
        call ga_copy_patch('n', 
1009
 
     $     g_r2, 1, n, 1, nvec, 
1010
 
     $     g_r(1), 1, n, 1,nvec)
1011
 
        if (ncomp.gt.1) call ga_copy_patch('n', 
1012
 
     $     g_r2, n+1, n2, 1, nvec, 
1013
 
     $     g_r(2), 1, n, 1,nvec)
1014
 
        if (lifetime) then
1015
 
          call ga_copy_patch('n', 
1016
 
     $       g_r2, n2+1, n3, 1, nvec, 
1017
 
     $       g_r_im(1), 1, n, 1,nvec)
1018
 
          if (ncomp.gt.1) call ga_copy_patch('n', 
1019
 
     $       g_r2, n3+1, n4, 1, nvec, 
1020
 
     $       g_r_im(2), 1, n, 1,nvec)
1021
 
        endif                   ! lifetime
1022
 
        call ga_sync()
1023
 
        do ipm = 1,ncomp
1024
 
          call ga_add(1.0d0, g_r(ipm), 1.0d0, g_x(ipm), g_x(ipm))
1025
 
          if (lifetime)
1026
 
     &       call ga_add(1.0d0, g_r_im(ipm),
1027
 
     &       1.0d0, g_x_im(ipm), g_x_im(ipm))
1028
 
        enddo
1029
 
        call ga_sync()
1030
 
        if (debug) write (6,*) 'lkain3_damp: gr to gx complete'
1031
 
c        
 
1311
       call update_g_x1(g_r2,    ! in   :
 
1312
     &                 g_x,     ! in/ou: updated
 
1313
     &                 g_x_im,  ! in/ou: updated
 
1314
     &                 ncomp,   ! in   : nr. components
 
1315
     &                 nvec,    ! in   : (x,y,z) = 3
 
1316
     &                 n,       ! in   : size of block (e.g. RE1,RE2,IM1 or IM2)
 
1317
     &                 lifetime)! in   : = .true. if complex
 
1318
 
 
1319
        if (debug1) then
 
1320
        if (ga_nodeid().eq.0)
 
1321
     &   write(*,*) '---------g_x-AFT-2(',iter,')-----START' 
 
1322
        call ga_print(g_x)
 
1323
        if (ga_nodeid().eq.0)
 
1324
     &   write(*,*) '---------g_x-AFT-2(',iter,')-----END'  
 
1325
        endif ! end-if-debug1
 
1326
       
1032
1327
        if (.not. ga_destroy(g_a)) call errquit
1033
1328
     &     ('lkain_2cpl: a',0, GA_ERR)
1034
1329
        if (.not. ga_destroy(g_bb)) call errquit
1038
1333
c       
1039
1334
c       Reduce the subspace as necessary
1040
1335
c       
 
1336
        if (ga_nodeid().eq.0) then
 
1337
         write(*,17) iter,nsub,maxsub,nvec,n4
 
1338
 17      format('(iter,nsub,maxsub,nvec,n4)=(',
 
1339
     &          i3,',',i4,',',i4,',',i3,',',i4,')')
 
1340
        endif
 
1341
        if (debug1) then
 
1342
         if (iter.eq.2) then
 
1343
          write(*,*) 'STOP-test-linsolver'
 
1344
          stop
 
1345
         endif
 
1346
        endif ! end-if-debug1 
 
1347
 
1041
1348
        if (nsub .eq. maxsub) then
1042
 
          do isub = nvec+1, maxsub, nvec
1043
 
c           Real part:
1044
 
c           component 1: 
1045
 
            call ga_copy_patch('n', 
1046
 
     $         g_Ay,    1, n, isub, isub+nvec-1, 
1047
 
     $         g_Ax(1), 1, n, 1, nvec)
1048
 
            call ga_sync()
1049
 
            call ga_copy_patch('n', 
1050
 
     $         g_Ax(1), 1, n, 1, nvec,
1051
 
     $         g_Ay,    1, n, isub-nvec, isub-1)
1052
 
c           component 2: 
1053
 
            if (ncomp.gt.1) then
1054
 
              call ga_copy_patch('n', 
1055
 
     $           g_Ay,    n+1, n2, isub, isub+nvec-1, 
1056
 
     $           g_Ax(2), 1,   n,  1, nvec)
1057
 
              call ga_sync()
1058
 
              call ga_copy_patch('n', 
1059
 
     $           g_Ax(2), 1,   n,  1, nvec,
1060
 
     $           g_Ay,    n+1, n2, isub-nvec, isub-1)  
1061
 
            endif
1062
 
            if (lifetime) then
1063
 
c             Imaginary part:
1064
 
c             component 1: 
1065
 
              call ga_copy_patch('n', 
1066
 
     $           g_Ay,    n2+1, n3, isub, isub+nvec-1, 
1067
 
     $           g_Ax_im(1), 1, n, 1, nvec)
1068
 
              call ga_sync()
1069
 
              call ga_copy_patch('n', 
1070
 
     $           g_Ax_im(1), 1, n, 1, nvec,
1071
 
     $           g_Ay,    n2+1, n3, isub-nvec, isub-1)
1072
 
c             component 2: 
1073
 
              if (ncomp.gt.1) then
1074
 
                call ga_copy_patch('n', 
1075
 
     $             g_Ay,    n3+1, n4, isub, isub+nvec-1, 
1076
 
     $             g_Ax_im(2), 1,   n,  1, nvec)
1077
 
                call ga_sync()
1078
 
                call ga_copy_patch('n', 
1079
 
     $             g_Ax_im(2), 1,   n,  1, nvec,
1080
 
     $             g_Ay,    n3+1, n4, isub-nvec, isub-1)         
1081
 
              endif
1082
 
            endif               ! lifetime
1083
 
c             
1084
 
c           Real part:
1085
 
c           component 1:
1086
 
            call ga_sync()
1087
 
            call ga_copy_patch('n', 
1088
 
     $         g_y,     1, n, isub, isub+nvec-1, 
1089
 
     $         g_Ax(1), 1, n, 1, nvec)
1090
 
            call ga_sync()
1091
 
            call ga_copy_patch('n', 
1092
 
     $         g_Ax(1), 1, n, 1, nvec,
1093
 
     $         g_y,     1, n, isub-nvec, isub-1)
1094
 
c           component 2:
1095
 
            if (ncomp.gt.1)  then
1096
 
              call ga_copy_patch('n', 
1097
 
     $           g_y,     n+1, n2, isub, isub+nvec-1, 
1098
 
     $           g_Ax(2), 1,   n,  1, nvec)
1099
 
              call ga_sync()
1100
 
              call ga_copy_patch('n', 
1101
 
     $           g_Ax(2), 1,   n,  1, nvec,
1102
 
     $           g_y,     n+1, n2, isub-nvec, isub-1)
1103
 
            endif
1104
 
c
1105
 
            if (lifetime) then
1106
 
c             Imaginary Part:
1107
 
c             component 1:
1108
 
              call ga_copy_patch('n', 
1109
 
     $           g_y,     n2+1, n3, isub, isub+nvec-1, 
1110
 
     $           g_Ax_im(1), 1, n, 1, nvec)
1111
 
              call ga_sync()
1112
 
              call ga_copy_patch('n', 
1113
 
     $           g_Ax_im(1), 1, n, 1, nvec,
1114
 
     $           g_y,     n2+1, n3, isub-nvec, isub-1)
1115
 
c             component 2:
1116
 
              if (ncomp.gt.1) then
1117
 
                call ga_copy_patch('n', 
1118
 
     $             g_y,     n3+1, n4, isub, isub+nvec-1, 
1119
 
     $             g_Ax_im(2), 1,   n,  1, nvec)
1120
 
                call ga_sync()
1121
 
                call ga_copy_patch('n', 
1122
 
     $             g_Ax_im(2), 1,   n,  1, nvec,
1123
 
     $             g_y,     n3+1, n4, isub-nvec, isub-1)
1124
 
              endif
1125
 
            endif               ! lifetime
1126
 
c           
1127
 
          end do                ! isub = nvec+1, maxsub, nvec
1128
 
          nsub = nsub - nvec
 
1349
c ====== FA: left-shifting patch ==== START
 
1350
c Note.- matrices Ay,y shift to left nvec positions
 
1351
c        removing leftmost patch of dimension: n4 x nvec
 
1352
         if (ga_nodeid().eq.0)
 
1353
     &    write(*,*) 'FA-matrix-nvec-left-shifting:'
 
1354
         do isub = nvec+1, maxsub, nvec
 
1355
          call ga_copy_patch('n',g_Ay,1,n4,isub,isub+nvec-1, 
 
1356
     $                           g_Ay,1,n4,isub-nvec,isub-1)
 
1357
          call ga_copy_patch('n',g_y ,1,n4,isub,isub+nvec-1, 
 
1358
     $                           g_y ,1,n4,isub-nvec,isub-1)
 
1359
c ====== FA: left-shifting patch ==== END
 
1360
         enddo ! end-loop-isub
 
1361
         nsub = nsub - nvec
1129
1362
        end if                  ! (nsub .eq. maxsub) 
1130
 
        if (debug) write (6,*) 'lkain3_damp: sub reduction comp.'
1131
 
c       
1132
1363
      enddo                     ! iter = 1,maxiter
1133
1364
  100 continue                  ! jump here if converged
1134
1365
c     
1151
1382
     &       ('lkain_2cpl: destroy',201, GA_ERR)
1152
1383
          if (.not. ga_destroy(g_r_im(ipm))) call errquit
1153
1384
     &       ('lkain_2cpl: destroy',51, GA_ERR)
1154
 
        endif
1155
 
c       
 
1385
        endif      
1156
1386
      enddo                     ! ipm = 1,2
1157
1387
      
1158
1388
      if (.not. ga_destroy(g_Ay)) call errquit
1170
1400
          call ga_print(g_x_im(1))
1171
1401
          if (ncomp.gt.1) call ga_print(g_x_im(2))
1172
1402
        endif
1173
 
      endif
1174
 
      
 
1403
      endif    
1175
1404
      
1176
1405
c ... jochen: disable this error exit during debuging phase
1177
1406
c     but print a warning instead