← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/util/ga_it2.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_it2.F
1
 
c $Id: ga_it2.F 19707 2010-10-29 17:59:36Z d3y133 $
 
1
c $Id: ga_it2.F 23368 2013-01-04 18:58:01Z niri $
 
2
C> \file ga_it2.F
 
3
C> \defgroup kain Krylov subspace Accelerated Inexact Newton method (KAIN)
 
4
C> \brief Implementation of Krylov-subspace nonlinear equation solver
 
5
C>
 
6
C> This files provides an implementation of the KAIN method to solve
 
7
C> linear and nonlinear systems of equations [1]. The solvers 
 
8
C> `ga_kain` and `ga_lkain` are general purpose routines. They require
 
9
C> a routine to calculate a matrix-vector product and a routine
 
10
C> to apply a preconditioner as arguments.
 
11
C>
 
12
C> [1] R.J. Harrison,
 
13
C>     <i>"Krylov Subspace Accelerated Inexact Newton Method for Linear
 
14
C>     and Nonlinear Systems",</i>
 
15
C>     Journal of Computational Chemistry (2004) <b>25</b>, pp 328-334,
 
16
C>     DOI: <a href="http://dx.doi.org/10.1002/jcc.10108">10.1002/jcc.10108</a>
 
17
C>
 
18
c
 
19
C> \ingroup kain
 
20
C> @{
 
21
C>
 
22
C> \brief A trivial matrix-vector product routine
 
23
C>
 
24
C> This routine is simply there to test the code for system solvers.
 
25
C> The matrix-vector product it calculates is simply a product involving
 
26
C> an explicitly stored matrix. This matrix is passed through a common
 
27
C> block.
 
28
C>
2
29
      subroutine test_product(acc,g_x, g_Ax)
3
30
      implicit none
4
 
      double precision acc
5
 
      integer g_x, g_Ax
 
31
#include "errquit.fh"
 
32
      double precision acc !< [Input] Accuracy (not used here)
 
33
      integer g_x  !< [Input] Global array with vector
 
34
      integer g_Ax !< [Output] Global array with matrix-vector product
6
35
c
7
36
      integer g_a, g_b
8
37
      common /testme/ g_A, g_B
9
 
      integer n, nvec, type
10
 
c
11
 
      call ga_inquire(g_x, type, n, nvec)
 
38
      integer n, nvec, typex, typeax
 
39
      double complex one, zero
 
40
c
 
41
      one  = cmplx(1.0d0,0.0d0)
 
42
      zero = cmplx(0.0d0,0.0d0)
 
43
c
 
44
      call ga_inquire(g_Ax, typeax, n, nvec)
 
45
      call ga_inquire(g_x, typex, n, nvec)
 
46
      if (typex.ne.typeax)
 
47
     $  call errquit("test_product: g_x not same type as g_Ax",
 
48
     $               typex,UERR)
12
49
      call ga_zero(g_Ax)
13
 
      call ga_dgemm('n', 'n', n, nvec, n, 1.0d0, g_A, g_x, 0.0d0, g_Ax)
 
50
      call ga_dgemm('n', 'n', n, nvec, n, one, g_A, g_x, zero, g_Ax)
14
51
c
15
52
      end
 
53
C>
 
54
C> \brief A trivial non-linear matrix-vector product routine
 
55
C>
 
56
C> This routine is simply there to test the code for system solvers.
 
57
C> The matrix-vector product it calculates is simply a product involving
 
58
C> an explicitly stored matrix and it adds another vector to achieve
 
59
C> something non-linear. The additional matrix and vector are passed
 
60
C> through a common block.
 
61
C>
16
62
      subroutine test_nlproduct(acc,g_x, g_Ax)
17
63
      implicit none
 
64
#include "errquit.fh"
18
65
      double precision acc
19
66
      integer g_x, g_Ax
20
67
c
21
68
      integer g_a, g_b
22
69
      common /testme/ g_A, g_B
23
 
      integer n, nvec, type
24
 
c
25
 
      call ga_inquire(g_x, type, n, nvec)
 
70
      integer n, nvec, typex, typeax, typeb
 
71
      double complex mone, one, zero
 
72
c
 
73
      mone = cmplx(-1.0d0,0.0d0)
 
74
      one  = cmplx(1.0d0,0.0d0)
 
75
      zero = cmplx(0.0d0,0.0d0)
 
76
c
 
77
      call ga_inquire(g_b, typeb, n, nvec)
 
78
      call ga_inquire(g_Ax, typeax, n, nvec)
 
79
      call ga_inquire(g_x, typex, n, nvec)
 
80
      if (typex.ne.typeax)
 
81
     $  call errquit("test_nlproduct: g_x not same type as g_Ax",
 
82
     $               typex,UERR)
 
83
      if (typex.ne.typeb)
 
84
     $  call errquit("test_nlproduct: g_x not same type as g_b",
 
85
     $               typex,UERR)
26
86
      call ga_zero(g_Ax)
27
 
      call ga_dgemm('n', 'n', n, nvec, n, 1.0d0, g_A, g_x, 0.0d0, g_Ax)
28
 
      call ga_add(1.0d0, g_AX, -1.0d0, g_B, g_AX)
 
87
      call ga_dgemm('n', 'n', n, nvec, n, one, g_A, g_x, zero, g_Ax)
 
88
      call ga_add(one, g_AX, mone, g_B, g_AX)
29
89
c
30
90
      end
 
91
C>
 
92
C> \brief A trivial preconditioner
 
93
C>
 
94
C> This routine is simply there to test the code for system solvers.
 
95
C> It simply scales the elements of a vector by
 
96
C> \f{eqnarray*}{
 
97
C>    x_i = \frac{x_i}{i - \epsilon}
 
98
C> \f}
 
99
C> where \f$ \epsilon \f$ is the shift.
 
100
C>
31
101
      subroutine test_precond(g_x,shift)
32
102
      implicit none
 
103
#include "mafdecls.fh"
 
104
#include "errquit.fh"
33
105
#include "global.fh"
34
106
      integer g_x
35
107
      double precision shift
36
108
      integer n, nvec, type
37
109
      integer i, ivec
38
110
      double precision x
 
111
      double complex zx
39
112
c
40
113
      call ga_inquire(g_x, type, n, nvec)
41
114
      if (ga_nodeid() .eq. 0) then
42
 
         do ivec = 1, nvec
43
 
            do i = 1, n
44
 
               call ga_get(g_x, i, i, ivec, ivec, x, 1)
45
 
               x = x / (dble(i) - shift)
46
 
               call ga_put(g_x, i, i, ivec, ivec, x, 1)
47
 
            end do
48
 
         end do
 
115
         if (type.eq.MT_DBL) then
 
116
            do ivec = 1, nvec
 
117
               do i = 1, n
 
118
                  call ga_get(g_x, i, i, ivec, ivec, x, 1)
 
119
                  x = x / (dble(i) - shift)
 
120
                  call ga_put(g_x, i, i, ivec, ivec, x, 1)
 
121
               end do
 
122
            end do
 
123
         else if (type.eq.MT_DCPL) then
 
124
            do ivec = 1, nvec
 
125
               do i = 1, n
 
126
                  call ga_get(g_x, i, i, ivec, ivec, zx, 1)
 
127
                  zx = zx / (dble(i) - shift)
 
128
                  call ga_put(g_x, i, i, ivec, ivec, zx, 1)
 
129
               end do
 
130
            end do
 
131
         else
 
132
            call errquit('test_precond: illegal type',type,UERR)
 
133
         endif
49
134
      end if
50
135
      call ga_sync()
51
136
c
52
137
      end
 
138
C>
 
139
C> \brief test driver for ga_lkain
 
140
C>
 
141
C> This routine drives a standard test of the ga_lkain 
 
142
C> and the ga_kain solver.
 
143
C> Because the test is fixed the subroutine takes no inputs
 
144
C> and it prints a summary on standard output. 
 
145
C>
53
146
      subroutine ga_lkain_test()
54
147
      implicit none
55
148
#include "errquit.fh"
56
149
#include "mafdecls.fh"
57
150
#include "global.fh"
58
151
      integer n, nvec, maxsub, maxiter
59
 
      parameter (n=300, nvec=1)
 
152
c     parameter (n=300, nvec=1)
 
153
      parameter (n=10, nvec=1)
 
154
      double complex angle
60
155
      integer g_x
61
 
      external test_precond, test_product, test_nlproduct
 
156
      double precision util_random
 
157
      external test_precond, test_product, test_nlproduct, util_random
62
158
c
 
159
      integer g_tx, g_ta, g_tb
63
160
      integer g_a, g_b
64
161
      common /testme/ g_a, g_b
65
162
c
66
163
      integer i
 
164
      logical  ga_copy_dz
 
165
      external ga_copy_dz
67
166
****      integer info
68
167
****      double precision a(n,n), b(n), w(n)
69
168
c
72
171
c
73
172
      if (.not. ga_create(MT_DBL, n, nvec, 'testx', 0, 0, g_x))
74
173
     $     call errquit('test kain', 1, GA_ERR)
75
 
      if (.not. ga_create(MT_DBL, n, nvec, 'testx', 0, 0, g_b))
 
174
      if (.not. ga_create(MT_DBL, n, nvec, 'testb', 0, 0, g_b))
76
175
     $     call errquit('test kain', 2, GA_ERR)
77
 
      if (.not. ga_create(MT_DBL, n, n, 'testx', 0, 0, g_A))
 
176
      if (.not. ga_create(MT_DBL, n, n, 'testA', 0, 0, g_A))
78
177
     $     call errquit('test kain', 3, GA_ERR)
79
178
c
80
179
      call ga_ran_fill(g_A, 1, n, 1, n)
88
187
c
89
188
      call ga_copy(g_b, g_x)
90
189
      call test_precond(g_x,0.0d0)
91
 
c
92
 
****      call ga_get(g_a, 1, n, 1, n, a, n)
93
 
****      call ga_get(g_b, 1, n, 1, nvec, b, n)
94
 
****      call dgesv(n, nvec, a, n, w, b, n, info)
95
 
****      write(6,*) ' info ', info
96
 
****      call ga_put(g_x, 1, n, 1, nvec, b, n)
97
 
c
98
 
C     This should have something other than zero
99
 
      call ga_lkain(0,g_x, g_b,test_product,test_precond,1d-6,maxsub,
100
 
     $     maxiter,.true.,.true.)
101
 
      call ga_summarize(0)
102
 
      call ma_summarize_allocated_blocks()
103
 
c
104
 
      write(6,*)
105
 
      write(6,*)
106
 
      write(6,*)
107
 
      write(6,*) ' DOING NON LINEAR '
108
 
      write(6,*)
109
 
      write(6,*)
110
 
      write(6,*)
111
 
      call ga_copy(g_b, g_x)
112
 
      call test_precond(g_x,0.0d0)
113
 
      maxsub = 10
114
 
c
115
 
      call ga_kain(g_x, 
116
 
     $     test_nlproduct, test_precond, 
117
 
     $     1d-6, 
118
 
     $     10.0d0, 10.0d0,
119
 
     $     maxsub, maxiter, 
120
 
     $     .true.)
121
 
c
 
190
      call util_flush(6)
 
191
      call ga_print(g_b)
 
192
      call ga_print(g_x)
 
193
      call util_flush(6)
 
194
c
 
195
****      call ga_get(g_a, 1, n, 1, n, a, n)
 
196
****      call ga_get(g_b, 1, n, 1, nvec, b, n)
 
197
****      call dgesv(n, nvec, a, n, w, b, n, info)
 
198
****      write(6,*) ' info ', info
 
199
****      call ga_put(g_x, 1, n, 1, nvec, b, n)
 
200
c
 
201
C     This should have something other than zero
 
202
      call ga_lkain(0,g_x, g_b,test_product,test_precond,1d-6,maxsub,
 
203
     $     maxiter,.true.,.true.)
 
204
      call util_flush(6)
 
205
      call ga_print(g_x)
 
206
      call util_flush(6)
 
207
      call ga_summarize(0)
 
208
      call ma_summarize_allocated_blocks()
 
209
c
 
210
      write(6,*)
 
211
      write(6,*)
 
212
      write(6,*)
 
213
      write(6,*) ' DOING NON LINEAR '
 
214
      write(6,*)
 
215
      write(6,*)
 
216
      write(6,*)
 
217
      call ga_copy(g_b, g_x)
 
218
      call test_precond(g_x,0.0d0)
 
219
      maxsub = 10
 
220
c
 
221
      call ga_kain(g_x, 
 
222
     $     test_nlproduct, test_precond, 
 
223
     $     1d-6, 
 
224
     $     10.0d0, 10.0d0,
 
225
     $     maxsub, maxiter, 
 
226
     $     .true.)
 
227
c
 
228
      call util_flush(6)
 
229
      call ga_print(g_x)
 
230
      call util_flush(6)
 
231
      call ga_summarize(0)
 
232
      call ma_summarize_allocated_blocks()
 
233
c
 
234
c     Same but now COMPLEX
 
235
c
 
236
      g_tx = g_x
 
237
      g_tb = g_b
 
238
      g_tA = g_A
 
239
      maxsub = 4*nvec
 
240
      maxiter = 100
 
241
      if (.not. ga_create(MT_DCPL, n, nvec, 'testx', 0, 0, g_x))
 
242
     $     call errquit('test kain', 1, GA_ERR)
 
243
      if (.not. ga_create(MT_DCPL, n, nvec, 'testb', 0, 0, g_b))
 
244
     $     call errquit('test kain', 2, GA_ERR)
 
245
      if (.not. ga_create(MT_DCPL, n, n, 'testA', 0, 0, g_A))
 
246
     $     call errquit('test kain', 3, GA_ERR)
 
247
c
 
248
c     call ga_ran_fill(g_A, 1, n, 1, n)
 
249
c     call ga_ran_fill(g_b, 1, n, 1, nvec)
 
250
c     if (ga_nodeid() .eq. 0) then
 
251
c        do i = 1, n
 
252
c           call ga_put(g_a, i, i, i, i, 0.5*cmplx(dble(i)), 1)
 
253
c        end do
 
254
c     end if
 
255
      if (.not. ga_copy_dz(g_tb,g_b)) 
 
256
     +  call errquit('test kain: copy g_b failed',0,GA_ERR)
 
257
      if (.not. ga_copy_dz(g_tA,g_A))
 
258
     +  call errquit('test kain: copy g_A failed',0,GA_ERR)
 
259
      call ga_sync()
 
260
      angle = cmplx(1.0d0,util_random(0))
 
261
      angle = angle/abs(angle)
 
262
      call ga_scale(g_b,angle)
 
263
      angle = cmplx(1.0d0,util_random(0))
 
264
      angle = angle/abs(angle)
 
265
      call ga_scale(g_A,angle)
 
266
c
 
267
      if (.not. ga_destroy(g_tx))
 
268
     $     call errquit('test kain', 1, GA_ERR)
 
269
      if (.not. ga_destroy(g_tb))
 
270
     $     call errquit('test kain', 2, GA_ERR)
 
271
      if (.not. ga_destroy(g_tA))
 
272
     $     call errquit('test kain', 3, GA_ERR)
 
273
c
 
274
      call ga_copy(g_b, g_x)
 
275
      call test_precond(g_x,0.0d0)
 
276
      call util_flush(6)
 
277
      call ga_print(g_b)
 
278
      call ga_print(g_x)
 
279
      call util_flush(6)
 
280
c
 
281
****      call ga_get(g_a, 1, n, 1, n, a, n)
 
282
****      call ga_get(g_b, 1, n, 1, nvec, b, n)
 
283
****      call dgesv(n, nvec, a, n, w, b, n, info)
 
284
****      write(6,*) ' info ', info
 
285
****      call ga_put(g_x, 1, n, 1, nvec, b, n)
 
286
c
 
287
C     This should have something other than zero
 
288
      call ga_lkain(0,g_x, g_b,test_product,test_precond,1d-6,maxsub,
 
289
     $     maxiter,.true.,.true.)
 
290
      call util_flush(6)
 
291
      call ga_print(g_x)
 
292
      call util_flush(6)
 
293
      call ga_summarize(0)
 
294
      call ma_summarize_allocated_blocks()
 
295
c
 
296
      write(6,*)
 
297
      write(6,*)
 
298
      write(6,*)
 
299
      write(6,*) ' DOING NON LINEAR '
 
300
      write(6,*)
 
301
      write(6,*)
 
302
      write(6,*)
 
303
      call ga_copy(g_b, g_x)
 
304
      call test_precond(g_x,0.0d0)
 
305
      maxsub = 10
 
306
c
 
307
      call ga_kain(g_x, 
 
308
     $     test_nlproduct, test_precond, 
 
309
     $     1d-6, 
 
310
     $     10.0d0, 10.0d0,
 
311
     $     maxsub, maxiter, 
 
312
     $     .true.)
 
313
c
 
314
      call util_flush(6)
 
315
      call ga_print(g_x)
 
316
      call util_flush(6)
122
317
      call ga_summarize(0)
123
318
      call ma_summarize_allocated_blocks()
124
319
      call errquit('done',0, MEM_ERR)
125
320
c
 
321
      if (.not. ga_destroy(g_x))
 
322
     $     call errquit('test kain', 1, GA_ERR)
 
323
      if (.not. ga_destroy(g_b))
 
324
     $     call errquit('test kain', 2, GA_ERR)
 
325
      if (.not. ga_destroy(g_A))
 
326
     $     call errquit('test kain', 3, GA_ERR)
 
327
c
126
328
      end
 
329
C>
 
330
C> \brief The Linear System Solver
 
331
C>
 
332
C> The routine solves a linear system of equations using a Krylov
 
333
C> subspace method:
 
334
C> \f{eqnarray*}{
 
335
C>   F(x) &=& b \\\\
 
336
C>   F(x) &=& Ax
 
337
C> \f}
 
338
C> The Right-Hand Side (RHS) matrix and the solution
 
339
C> matrix store the data pertaining to a single linear system in 
 
340
C> columns. Using multiple columns the solver can solve multiple
 
341
C> equations simultaneously. Internally one Krylov subspace is used
 
342
C> for all systems of equations, and subspace vectors contribute to
 
343
C> the solutions of all equations.
 
344
C>
 
345
C> An special argument `odiff` is provided that allows one to use a
 
346
C> difference strategy for the required matrix-vector products. I.e.
 
347
C> rather than computing 
 
348
C> \f{eqnarray*}{
 
349
C>     y &=& Ax_i
 
350
C> \f}
 
351
C> each iteration one uses a formulation that requires
 
352
C> \f{eqnarray*}{
 
353
C>     y' &=& A(x_i-x_{i-1})
 
354
C> \f}
 
355
C> As \f$ x \f$ converges \f$ x_i - x_{i-1} \f$ approaches \f$ 0 \f$.
 
356
C> This may be exploited in the matrix-vector product to achieve 
 
357
C> considerable savings.
 
358
C>
 
359
C> Furthermore the routine takes two subroutines as arguments. One 
 
360
C> provides the matrix-vector product and has the interface
 
361
C> ~~~~
 
362
C> subroutine product(acc,g_x,g_Ax)
 
363
C> ~~~~
 
364
C> which evaluates
 
365
C> \f{eqnarray*}{
 
366
C>   g_{Ax} &=& A g_x
 
367
C> \f}
 
368
C> where
 
369
C>
 
370
C> - acc: the accuracy required for the matrix-vector elements
 
371
C>
 
372
C> - g_x: is a global array containing the input vectors in columns
 
373
C> 
 
374
C> - g_Ax: is a global array that will contain the matrix-vector product
 
375
C>   results in columns
 
376
C>
 
377
C> Note that `ga_lkain` may pass global arrays with a different number
 
378
C> of vectors than the number of linear systems being solved. Hence
 
379
C> the product routine should check the actual dimension of g_x and
 
380
C> g_Ax.
 
381
C>
 
382
C> The second subroutine provides a preconditioner and has the 
 
383
C> interface
 
384
C> ~~~~
 
385
C> subroutine precond(g_x,shift)
 
386
C> ~~~~
 
387
C> where 
 
388
C> 
 
389
C> - g_x: is a global array holding a set of vectors that will have
 
390
C>   the preconditioner applied to them 
 
391
C>
 
392
C> - shift: is a shift to ensure the preconditioner is non-singular
 
393
C>
 
394
C> Note that the number of vectors may also here be different from the
 
395
C> number of linear systems being solved. Hence the dimension of g_x
 
396
C> should be tested. Finally, the shift is not used in this linear
 
397
C> system solver. The argument is provided to ensure that the same
 
398
C> preconditioner can be used in eigenvalue solvers where the shift is
 
399
C> used.
 
400
C> 
127
401
      subroutine ga_lkain(rtdb,g_x, g_b, product, precond, 
128
402
     $     tol, mmaxsub, maxiter, odiff, oprint)
129
403
      implicit none
133
407
#include "util.fh"
134
408
#include "rtdb.fh"
135
409
c
136
 
      integer g_x               ! [input/output] Initial guess/solution
137
 
      integer g_b               ! [input] Right-hand side vectors
138
 
      external product          ! [input] product routine
139
 
      external precond          ! [input] preconditioner routine
140
 
      double precision tol      ! [input] convergence threshold
141
 
      integer mmaxsub           ! [input] maximum subspace dimension
142
 
      integer maxiter           ! [input] maximum no. of iterations
143
 
      logical odiff             ! [input] use differences in product
144
 
      logical oprint            ! [input] print flag
145
 
      integer rtdb
 
410
      integer g_x          !< [Input/Output] Initial guess/solution
 
411
      integer g_b          !< [Input] Right-hand side vectors
 
412
      external product     !< [Input] product routine
 
413
      external precond     !< [Input] preconditioner routine
 
414
      double precision tol !< [Input] convergence threshold
 
415
      integer mmaxsub      !< [Input] maximum subspace dimension
 
416
      integer maxiter      !< [Input] maximum no. of iterations
 
417
      logical odiff        !< [Input] use differences in product
 
418
      logical oprint       !< [Input] print flag
 
419
      integer rtdb         !< [Input] the RTDB handle
146
420
c
147
 
c     Solves the linear equations A(X)=0 for multiple vectors.
 
421
c     Solves the linear equations A(X)-B=0 for multiple vectors.
148
422
c
149
423
c     call product(acc,g_x, g_Ax)
150
424
c     . acc is the accuracy trequired for each element of the product
165
439
c     Needs to be extended to store the sub-space vectors out-of-core
166
440
c     at least while the product() routine is being executed.
167
441
c
168
 
      integer iter, n, nvec, nsub, isub, type, maxsub
 
442
      integer dimi, dimj
 
443
      integer iter, n, nvec, nsub, isub, typex, typeb, maxsub
169
444
      integer g_y, g_Ay, g_Ax, g_r, g_a, g_bb, g_c, g_xold, g_Axold
170
445
      double precision rmax, acc, ga_svd_tol
 
446
      double complex zero, one, mone
171
447
      logical converged
172
448
      logical odebug
173
449
c
 
450
      zero = cmplx(0.0d0,0.0d0)
 
451
      one  = cmplx(1.0d0,0.0d0)
 
452
      mone = cmplx(-1.0d0,0.0d0)
 
453
c
174
454
      odebug = util_print('debug lsolve', print_never) .and. 
175
455
     $     ga_nodeid().eq.0
176
456
      if (.not.rtdb_get(rtdb,'cphf:acc',mt_dbl,1,acc)) acc=0.01d0*tol
177
 
      call ga_inquire(g_x, type, n, nvec)
 
457
      call ga_inquire(g_b, typeb, dimi, dimj)
 
458
      call ga_inquire(g_x, typex, n, nvec)
 
459
      if (typeb.ne.typex) 
 
460
     $  call errquit("ga_lkain: solution and right-hand-side vectors not
 
461
     $ of same type",0,UERR)
178
462
      maxsub = mmaxsub          ! So don't modify input scalar arg
179
463
      if (maxsub .lt. 3*nvec) maxsub = 3*nvec
180
464
      maxsub = (maxsub/nvec)*nvec
202
486
         call util_flush(6)
203
487
      end if
204
488
c
205
 
      if (.not. ga_create(MT_DBL, n, maxsub, 'lkain: Y', 
 
489
      if (.not. ga_create(typex, n, maxsub, 'lkain: Y', 
206
490
     $     0, 0, g_y))
207
491
     $     call errquit('lkain: failed allocating subspace', maxsub,
208
492
     &       GA_ERR)
209
 
      if (.not. ga_create(MT_DBL, n, maxsub, 'lkain: Ay', 
 
493
      if (.not. ga_create(typex, n, maxsub, 'lkain: Ay', 
210
494
     $     0, 0, g_Ay))
211
495
     $     call errquit('lkain: failed allocating subspace2', maxsub,
212
496
     &       GA_ERR)
213
 
      if (.not. ga_create(MT_DBL, n, nvec, 'lkain: Ax',
 
497
      if (.not. ga_create(typex, n, nvec, 'lkain: Ax',
214
498
     $     0, 0, g_Ax))
215
499
     $     call errquit('lkain: failed allocating subspace3', nvec,
216
500
     &       GA_ERR)
217
 
      if (.not. ga_create(MT_DBL, n, nvec, 'lkain: r',
 
501
      if (.not. ga_create(typex, n, nvec, 'lkain: r',
218
502
     $     0, 0, g_r))
219
503
     $     call errquit('lkain: failed allocating subspace4', nvec,
220
504
     &       GA_ERR)
221
505
      if (odiff) then
222
 
         if (.not. ga_create(MT_DBL, n, nvec, 'lkain: xold',
 
506
         if (.not. ga_create(typex, n, nvec, 'lkain: xold',
223
507
     $        0, 0, g_xold))
224
508
     $        call errquit('lkain: failed allocating subspace5', nvec,
225
509
     &       GA_ERR)
226
 
         if (.not. ga_create(MT_DBL, n, nvec, 'lkain: xold',
 
510
         if (.not. ga_create(typex, n, nvec, 'lkain: Axold',
227
511
     $        0, 0, g_Axold))
228
512
     $        call errquit('lkain: failed allocating subspace6', nvec,
229
513
     &       GA_ERR)
247
531
      converged = .false.
248
532
      do iter = 1, maxiter
249
533
         if (odiff) then
250
 
            call ga_add(1.0d0, g_x, -1.0d0, g_xold,  g_x)
 
534
            call ga_add(one, g_x, mone, g_xold,  g_x)
251
535
         end if
252
536
         call product(acc,g_x, g_Ax)
253
537
         if (odiff) then
254
 
            call ga_add(1.0d0, g_Ax, 1.0d0, g_Axold, g_Ax)
255
 
            call ga_add(1.0d0, g_x,  1.0d0, g_xold,  g_x)
 
538
            call ga_add(one, g_Ax, one, g_Axold, g_Ax)
 
539
            call ga_add(one, g_x,  one, g_xold,  g_x)
256
540
            call ga_copy(g_x, g_xold)
257
541
            call ga_copy(g_Ax, g_Axold)
258
542
         end if
259
543
         call ga_zero(g_r)
260
544
         call ga_sync()
261
 
         call ga_add(1.0d0, g_b, -1.0d0, g_Ax, g_r) ! The residual
 
545
         call ga_add(one, g_b, mone, g_Ax, g_r) ! The residual
262
546
         call ga_sync()
263
547
         call ga_maxelt(g_r, rmax)
264
548
         if (oprint .and. ga_nodeid().eq.0) then
287
571
c     Form and solve the subspace equations using SVD in order
288
572
c     to manage near linear dependence in the subspace.
289
573
c     
290
 
         if (.not. ga_create(MT_DBL, nsub, nsub, 'lkain: A', 0, 0, g_a))
 
574
         if (.not. ga_create(typex, nsub, nsub, 'lkain: A', 0, 0, g_a))
291
575
     $        call errquit('lkain: allocating g_a?', nsub, GA_ERR)
292
 
         if (.not. ga_create(MT_DBL, nsub, nvec, 'lkain: B', 0, 0,g_bb))
 
576
         if (.not. ga_create(typex, nsub, nvec, 'lkain: B', 0, 0,g_bb))
293
577
     $        call errquit('lkain: allocating g_bb?', nsub, GA_ERR)
294
 
         if (.not. ga_create(MT_DBL, nsub, nvec, 'lkain: C', 0, 0, g_c))
 
578
         if (.not. ga_create(typex, nsub, nvec, 'lkain: C', 0, 0, g_c))
295
579
     $        call errquit('lkain: allocating g_c?', nsub, GA_ERR)
296
580
         call ga_zero(g_a)
297
581
         call ga_zero(g_bb)
298
582
         call ga_zero(g_c)
299
583
         call ga_sync()
300
 
         call ga_dgemm('t','n',nsub,nsub,n,1.0d0,g_y,g_Ay,0.0d0,g_a)
301
 
         call ga_sync()
302
 
         call ga_dgemm('t','n',nsub,nvec,n,1.0d0,g_y,g_r,0.0d0,g_bb)
303
 
         call ga_sync()
304
 
         if (odebug) call ga_print(g_a)
305
 
         if (odebug) call ga_print(g_c)
 
584
         call ga_dgemm('c','n',nsub,nsub,n,one,g_y,g_Ay,zero,g_a)
 
585
         call ga_sync()
 
586
         call ga_dgemm('c','n',nsub,nvec,n,one,g_y,g_r,zero,g_bb)
 
587
         call ga_sync()
 
588
         if (odebug) then
 
589
           call util_flush(6)
 
590
           call ga_print(g_a)
 
591
           call ga_print(g_bb)
 
592
           call util_flush(6)
 
593
         endif
306
594
c
307
595
c     The threshold used here should reflect the accuracy in the
308
596
c     products.  If very accurate products are used, then there is big
313
601
c     screening with a realistic threshold is important.
314
602
c
315
603
         call ga_svd_solve_seq(g_a,g_bb,g_c,ga_svd_tol)
316
 
         if (odebug) call ga_print(g_c)
 
604
         if (odebug) then
 
605
           call util_flush(6)
 
606
           call ga_print(g_c)
 
607
           call util_flush(6)
 
608
         endif
317
609
c
318
610
c     Form and add the correction, in parts, onto the solution
319
611
c
320
612
         call ga_sync()
321
 
         call ga_dgemm('n','n',n,nvec,nsub,-1.0d0,g_Ay,g_c,1.0d0,g_r)
 
613
         call ga_dgemm('n','n',n,nvec,nsub,mone,g_Ay,g_c,one,g_r)
322
614
         if (odebug) then
 
615
            call util_flush(6)
323
616
            write(6,*) ' The update in the complement '
 
617
            call util_flush(6)
324
618
            call ga_print(g_r)
325
619
         end if
326
620
         call ga_sync()
327
 
         call ga_add(1.0d0, g_r, 1.0d0, g_x, g_x)
 
621
         call ga_add(one, g_r, one, g_x, g_x)
328
622
         call ga_sync()
329
 
         call ga_dgemm('n','n',n,nvec,nsub,1.0d0,g_y,g_c,0.0d0,g_r)
 
623
         call ga_dgemm('n','n',n,nvec,nsub,one,g_y,g_c,zero,g_r)
330
624
         if (odebug) then
 
625
            call util_flush(6)
331
626
            write(6,*) ' The update in the subspace '
 
627
            call util_flush(6)
332
628
            call ga_print(g_r)
333
629
         end if
334
630
         call ga_sync()
335
 
         call ga_add(1.0d0, g_r, 1.0d0, g_x, g_x)
 
631
         call ga_add(one, g_r, one, g_x, g_x)
336
632
         call ga_sync()
337
633
c
338
634
         if (.not. ga_destroy(g_a)) call errquit('lkain: a',0, GA_ERR)
383
679
     &       CALC_ERR)
384
680
c
385
681
      end
 
682
C>
 
683
C> \brief The Non-Linear System Solver
 
684
C>
 
685
C> The routine solves a non-linear system of equations using a Krylov
 
686
C> subspace method.
 
687
C> \f{eqnarray*}{
 
688
C>   F(x) &=& 0 
 
689
C> \f}
 
690
C> For linear systems the vector function in the above equation would
 
691
C> be \f$ F(x) = Ax - b \f$ but this routine also allows for non-linear
 
692
C> vector functions.
 
693
C> The Right-Hand Side (RHS) matrix and the solution
 
694
C> matrix store the data pertaining to a single linear system in 
 
695
C> columns. Using multiple columns the solver can solve multiple
 
696
C> equations simultaneously. Internally one Krylov subspace is used
 
697
C> for all systems of equations, and subspace vectors contribute to
 
698
C> the solutions of all equations.
 
699
C>
 
700
C> Furthermore the routine takes two subroutines as arguments. One 
 
701
C> provides the vector residual and has the interface
 
702
C> ~~~~
 
703
C> subroutine residual(acc,g_x,g_Ax)
 
704
C> ~~~~
 
705
C> where
 
706
C>
 
707
C> - acc: the accuracy required for the matrix-vector elements
 
708
C>
 
709
C> - g_x: is a global array containing the input vectors in columns
 
710
C> 
 
711
C> - g_Ax: is a global array that will contain the residual-vector
 
712
C>   results in columns
 
713
C>
 
714
C> Note that `ga_lkain` may pass global arrays with a different number
 
715
C> of vectors than the number of linear systems being solved. Hence
 
716
C> the product routine should check the actual dimension of g_x and
 
717
C> g_Ax.
 
718
C>
 
719
C> The second subroutine provides a preconditioner and has the 
 
720
C> interface
 
721
C> ~~~~
 
722
C> subroutine precond(g_x,shift)
 
723
C> ~~~~
 
724
C> where 
 
725
C> 
 
726
C> - g_x: is a global array holding a set of vectors that will have
 
727
C>   the preconditioner applied to them 
 
728
C>
 
729
C> - shift: is a shift to ensure the preconditioner is non-singular
 
730
C>
 
731
C> Note that the number of vectors may also here be different from the
 
732
C> number of linear systems being solved. Hence the dimension of g_x
 
733
C> should be tested. Finally, the shift is not used in this linear
 
734
C> system solver. The argument is provided to ensure that the same
 
735
C> preconditioner can be used in eigenvalue solvers where the shift is
 
736
C> used.
 
737
C> 
386
738
      subroutine ga_kain(
387
739
     $     g_x, 
388
 
     $     product, precond, 
 
740
     $     residual, precond, 
389
741
     $     tol, 
390
742
     $     trustmin, trustmax, 
391
743
     $     maxsub, maxiter, 
397
749
#include "tcgmsg.fh"
398
750
#include "util.fh"
399
751
c
400
 
      integer g_x               ! [input/output] Initial guess/solution
401
 
      external product          ! [input] product routine
402
 
      external precond          ! [input] preconditioner routine
403
 
      double precision tol      ! [input] convergence threshold
404
 
      double precision trustmin, trustmax ! [input] range to constrain trust radius
405
 
      integer maxsub           ! [input] maximum subspace dimension
406
 
      integer maxiter           ! [input] maximum no. of iterations
407
 
      logical oprint            ! [input] print flag
 
752
      integer g_x               !< [Input/Output] Initial guess/solution
 
753
      external residual         !< [Input] residual routine
 
754
      external precond          !< [Input] preconditioner routine
 
755
      double precision tol      !< [Input] convergence threshold
 
756
      double precision trustmin !< [Input] range to constrain trust radius
 
757
      double precision trustmax !< [Input] range to constrain trust radius
 
758
      integer maxsub            !< [Input] maximum subspace dimension
 
759
      integer maxiter           !< [Input] maximum no. of iterations
 
760
      logical oprint            !< [Input] print flag
408
761
c
409
762
c     Solves the non-linear equations f(X)=0 for multiple vectors.
410
763
c
411
 
c     call product(acc, g_x, g_Ax)
412
 
c     . acc is the accuracy trequired for each element of the product
 
764
c     call residual(acc, g_x, g_Ax)
 
765
c     . acc is the accuracy trequired for each element of the residual
413
766
c     . g_x contains the vectors and g_Ax should be filled
414
 
c     .     with the product vectors.  The no. of vectors (columns) in
 
767
c     .     with the residual vectors.  The no. of vectors (columns) in
415
768
c     . g_x might differ from the no. of vectors input to ga_kain().
416
769
c
417
770
c     call precond(g_x,shift)
423
776
c
424
777
      integer maxmaxsub
425
778
      parameter (maxmaxsub = 20)
426
 
      integer iter, n, nvec, nsub, isub, jsub, type
 
779
      integer iter, n, nvec, nsub, isub, jsub, typex
427
780
      integer g_y, g_Ay, g_Ax, g_delta, g_a, g_b, g_c
428
781
      double precision rmax, acc, trust
429
782
      double precision a(maxmaxsub,maxmaxsub), b(maxmaxsub), 
430
 
     $     c(maxmaxsub), csum
 
783
     $                 c(maxmaxsub), csum
 
784
      double complex za(maxmaxsub,maxmaxsub), zb(maxmaxsub), 
 
785
     $               zc(maxmaxsub), zcsum
 
786
      double complex zero, one, mone
431
787
      logical converged
432
788
      logical odebug
433
789
c
 
790
      zero = cmplx(0.0d0,0.0d0)
 
791
      one  = cmplx(1.0d0,0.0d0)
 
792
      mone = cmplx(-1.0d0,0.0d0)
 
793
c
434
794
      trust = trustmin
435
795
      acc = 0.01d0*tol
436
796
      if (maxsub .gt. maxmaxsub) maxsub = maxmaxsub
437
797
      odebug = util_print('debug lsolve', print_never) .and. 
438
798
     $     ga_nodeid().eq.0
439
799
c
440
 
      call ga_inquire(g_x, type, n, nvec)
 
800
      call ga_inquire(g_x, typex, n, nvec)
441
801
      if (nvec .ne. 1) call errquit('kain: nvec?', nvec, GA_ERR)
442
802
c
443
803
      if (oprint .and. ga_nodeid().eq.0) then
451
811
         call util_flush(6)
452
812
      end if
453
813
c
454
 
      if (.not. ga_create(MT_DBL, n, maxsub, 'kain: Y', 
 
814
      if (.not. ga_create(typex, n, maxsub, 'kain: Y', 
455
815
     $     0, 0, g_y))
456
816
     $     call errquit('kain: failed allocating subspace', maxsub,
457
817
     &       GA_ERR)
458
 
      if (.not. ga_create(MT_DBL, n, maxsub, 'kain: Ay', 
 
818
      if (.not. ga_create(typex, n, maxsub, 'kain: Ay', 
459
819
     $     0, 0, g_Ay))
460
820
     $     call errquit('kain: failed allocating subspace2', maxsub,
461
821
     &       GA_ERR)
462
 
      if (.not. ga_create(MT_DBL, n, 1, 'kain: Ax',
 
822
      if (.not. ga_create(typex, n, 1, 'kain: Ax',
463
823
     $     0, 0, g_Ax))
464
824
     $     call errquit('kain: failed allocating subspace3', 1, GA_ERR)
465
825
      call ga_zero(g_y)
475
835
      converged = .false.
476
836
      do iter = 1, maxiter
477
837
         call ga_zero(g_Ax)
478
 
         call product(acc, g_x, g_Ax)
 
838
         call residual(acc, g_x, g_Ax)
479
839
         call ga_maxelt(g_Ax, rmax)
480
840
         if (oprint .and. ga_nodeid().eq.0) then
481
841
            write(6,3) iter, nsub+1, rmax, util_wallsec()
500
860
c     Not converged ... make the update
501
861
c
502
862
         g_delta = g_Ax         ! A reminder that these two are aliased
503
 
         call ga_scale(g_delta, -1.0d0)
 
863
         call ga_scale(g_delta, mone)
504
864
c
505
865
         if (iter .gt. 1) then
506
866
c
507
867
c     Form the reduced space matrix and RHS
508
868
c
509
 
            call ga_local_mdot(n, nsub, nsub, a, maxmaxsub, g_y, g_Ay)
510
 
            do isub = 1, nsub-1
511
 
               b(isub) = -(a(isub,nsub) - a(nsub,nsub))
512
 
            end do
513
 
            do isub = 1, nsub-1
514
 
               do jsub = 1, nsub-1
515
 
                  a(isub,jsub) = a(isub,jsub)
516
 
     $                 - a(nsub,jsub) - a(isub,nsub) + a(nsub,nsub)
517
 
               end do
518
 
            end do
 
869
           if (typex.eq.MT_DBL) then
 
870
             call ga_local_mdot(n,nsub,nsub,a,maxmaxsub,g_y,g_Ay)
 
871
             do isub = 1, nsub-1
 
872
                b(isub) = -(a(isub,nsub) - a(nsub,nsub))
 
873
             end do
 
874
             do isub = 1, nsub-1
 
875
                do jsub = 1, nsub-1
 
876
                   a(isub,jsub) = a(isub,jsub)
 
877
     $                  - a(nsub,jsub) - a(isub,nsub) + a(nsub,nsub)
 
878
                end do
 
879
             end do
 
880
           else if (typex.eq.MT_DCPL) then
 
881
             call ga_local_zmdot(n,nsub,nsub,za,maxmaxsub,g_y,g_Ay)
 
882
             do isub = 1, nsub-1
 
883
                zb(isub) = -(za(isub,nsub) - za(nsub,nsub))
 
884
             end do
 
885
             do isub = 1, nsub-1
 
886
                do jsub = 1, nsub-1
 
887
                   za(isub,jsub) = za(isub,jsub)
 
888
     $                  - za(nsub,jsub) - za(isub,nsub) + za(nsub,nsub)
 
889
                end do
 
890
             end do
 
891
           else
 
892
             call errquit("ga_kain: illegal type",typex,UERR)
 
893
           endif
519
894
c
520
895
c     Solve the subspace equations (lazily using existing GA routine)
521
896
c     
522
 
            if (.not. ga_create(MT_DBL,nsub-1,nsub-1,'kain: A',
 
897
            if (.not. ga_create(typex,nsub-1,nsub-1,'kain: A',
523
898
     $           nsub-1,nsub-1,g_a))
524
899
     $           call errquit('kain: allocating g_a?', nsub, GA_ERR)
525
 
            if (.not. ga_create(MT_DBL,nsub-1,1,'kain: B',nsub-1,1,g_b))
 
900
            if (.not. ga_create(typex,nsub-1,1,'kain: B',nsub-1,1,g_b))
526
901
     $           call errquit('kain: allocating g_bb?', nsub, GA_ERR)
527
 
            if (.not. ga_create(MT_DBL,nsub-1,1,'kain: C',nsub-1,1,g_c))
 
902
            if (.not. ga_create(typex,nsub-1,1,'kain: C',nsub-1,1,g_c))
528
903
     $           call errquit('kain: allocating g_c?', nsub, GA_ERR)
529
904
            if (ga_nodeid() .eq. 0) then
530
 
               call ga_put(g_a, 1, nsub-1, 1, nsub-1, a, maxmaxsub)
531
 
               call ga_put(g_b, 1, nsub-1, 1, 1, b, 1)
 
905
              if (typex.eq.MT_DBL) then
 
906
                call ga_put(g_a, 1, nsub-1, 1, nsub-1, a, maxmaxsub)
 
907
                call ga_put(g_b, 1, nsub-1, 1, 1, b, 1)
 
908
              else if (typex.eq.MT_DCPL) then
 
909
                call ga_put(g_a, 1, nsub-1, 1, nsub-1, za, maxmaxsub)
 
910
                call ga_put(g_b, 1, nsub-1, 1, 1, zb, 1)
 
911
              else
 
912
                call errquit("ga_kain: illegal type",typex,UERR)
 
913
              endif
532
914
            end if
533
915
            call ga_sync
534
916
c
535
917
            call ga_svd_solve_seq(g_a,g_b,g_c,1d-14)
536
918
c
537
 
            if (odebug) call ga_print(g_c)
538
 
            if (ga_nodeid() .eq. 0) 
 
919
            if (odebug) then
 
920
              call util_flush(6)
 
921
              call ga_print(g_c)
 
922
              call util_flush(6)
 
923
            endif
 
924
            if (typex.eq.MT_DBL) then
 
925
              if (ga_nodeid() .eq. 0)
539
926
     $           call ga_get(g_c, 1, nsub-1, 1, 1, c, 1)
540
 
            call ga_brdcst(1, c, mdtob(nsub-1), 0)
541
 
            write(6,*) ' KAIN SUBSPACE COEFFS'
542
 
            call output(c, 1, nsub-1, 1, 1, nsub-1, 1, 1)
 
927
              call ga_brdcst(1, c, mdtob(nsub-1), 0)
 
928
              write(6,*) ' KAIN SUBSPACE COEFFS'
 
929
              call output(c, 1, nsub-1, 1, 1, nsub-1, 1, 1)
 
930
            else if (typex.eq.MT_DCPL) then
 
931
              if (ga_nodeid() .eq. 0)
 
932
     $           call ga_get(g_c, 1, nsub-1, 1, 1, zc, 1)
 
933
              call ga_brdcst(1, zc, 2*mdtob(nsub-1), 0)
 
934
              write(6,*) ' KAIN SUBSPACE COEFFS'
 
935
              call zoutput(zc, 1, nsub-1, 1, 1, nsub-1, 1, 1)
 
936
            endif
543
937
            call ga_sync
544
938
            if (.not. ga_destroy(g_a)) call errquit('kain: a',0, GA_ERR)
545
939
            if (.not. ga_destroy(g_b)) call errquit('kain: b',0, GA_ERR)
547
941
c     
548
942
c     Form the correction
549
943
c     
550
 
            csum = 0.0d0
551
 
            do isub = 1, nsub-1
552
 
               csum = csum + c(isub)
553
 
               call ga_add_patch( c(isub),  g_y, 1, n, isub, isub, 
554
 
     $              1.0d0, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
555
 
               call ga_add_patch(-c(isub), g_Ay, 1, n, isub, isub, 
556
 
     $              1.0d0, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
557
 
            end do
558
 
            call ga_add_patch(-csum,  g_y, 1, n, nsub, nsub, 
559
 
     $           1.0d0, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
560
 
            call ga_add_patch( csum,  g_Ay, 1, n, nsub, nsub, 
561
 
     $           1.0d0, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
944
            if (typex.eq.MT_DBL) then
 
945
              csum = 0.0d0
 
946
              do isub = 1, nsub-1
 
947
                 csum = csum + c(isub)
 
948
                 call ga_add_patch( c(isub),  g_y, 1, n, isub, isub, 
 
949
     $                one, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
950
                 call ga_add_patch(-c(isub), g_Ay, 1, n, isub, isub, 
 
951
     $                one, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
952
              end do
 
953
              call ga_add_patch(-csum,  g_y, 1, n, nsub, nsub, 
 
954
     $             one, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
955
              call ga_add_patch( csum,  g_Ay, 1, n, nsub, nsub, 
 
956
     $             one, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
957
            else if (typex.eq.MT_DCPL) then
 
958
              zcsum = 0.0d0
 
959
              do isub = 1, nsub-1
 
960
                 zcsum = zcsum + zc(isub)
 
961
                 call ga_add_patch( zc(isub),  g_y, 1, n, isub, isub, 
 
962
     $                one, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
963
                 call ga_add_patch(-zc(isub), g_Ay, 1, n, isub, isub, 
 
964
     $                one, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
965
              end do
 
966
              call ga_add_patch(-zcsum,  g_y, 1, n, nsub, nsub, 
 
967
     $             one, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
968
              call ga_add_patch( zcsum,  g_Ay, 1, n, nsub, nsub, 
 
969
     $             one, g_delta, 1, n, 1, 1, g_delta, 1, n, 1, 1)
 
970
            endif
562
971
         endif
563
972
c
564
973
c     Step restriction
569
978
            call ga_scale(g_delta, trust/rmax)
570
979
         end if
571
980
c
572
 
         call ga_add(1.0d0, g_delta, 1.0d0, g_x, g_x)
 
981
         call ga_add(one, g_delta, one, g_x, g_x)
573
982
c     
574
983
c     Reduce the subspace as necessary (note g_delta=g_Ax destroyed)
575
984
c     
605
1014
     &       CALC_ERR)
606
1015
c
607
1016
      end
 
1017
C>
 
1018
C> \brief A direct linear system solver
 
1019
C>
 
1020
C> Solve for X from the linear equations
 
1021
C> \f{eqnarray*}{
 
1022
C>    A*X &=& B
 
1023
C> \f}
 
1024
C> or more explicitly
 
1025
C> \f{eqnarray*}{
 
1026
C>    A(m,n)*X(n,nvec) = B(m,nvec)
 
1027
C> \f}
 
1028
C> Where \f$ A \f$ is a general real matrix (not necessarily square, or
 
1029
C> symmetric, or full rank) and \f$ X \f$ and \f$ B \f$ are matrices with one or more
 
1030
C> columns representing the solutions and right hand sides.  Singular
 
1031
C> values of \f$ A \f$ less than \f$ tol\f$  are neglected. \f$ X \f$ is returned.
 
1032
C>
 
1033
C> If the SVD of \f$ A \f$ is \f$ U*values*VT \f$, then the solution
 
1034
C> is of the form
 
1035
C> \f{eqnarray*}{
 
1036
C>    V*(1/values)*UT*B
 
1037
C> \f}
 
1038
C> where the reciprocal of `values` less than `tol` are neglected.
 
1039
C     
608
1040
      subroutine ga_svd_solve_seq(g_a, g_b, g_x, tol)
609
1041
      implicit none
610
1042
#include "errquit.fh"
611
1043
#include "global.fh"
612
1044
#include "mafdecls.fh"
613
1045
#include "util.fh"
614
 
      integer g_a, g_b, g_x
615
 
      double precision tol
 
1046
      integer g_a !< [Input] the matrix stored explicitly
 
1047
      integer g_b !< [Input] the right-hand-sides
 
1048
      integer g_x !< [Input/Output] the guess/solution
 
1049
      double precision tol !< [Input] the tolerance
616
1050
c
617
1051
c     Solve for X from the linear equations
618
1052
c
633
1067
c     where the reciprocal of values less than tol are neglected.
634
1068
c     
635
1069
      integer m,n,nn,type,nvec,nsing,l_val, k_val,g_u,g_vt,i,g_tmp
 
1070
      double complex one, zero
636
1071
      logical oprint
637
1072
c
638
1073
      oprint = util_print('debug svdsolve', print_high) .and.
639
1074
     $     ga_nodeid().eq.0
640
1075
c
 
1076
      one  = cmplx(1.0d0,0.0d0)
 
1077
      zero = cmplx(0.0d0,0.0d0)
 
1078
c
641
1079
      call ga_inquire(g_a, type, m, n)
642
1080
      call ga_inquire(g_b, type, nn, nvec)
643
1081
      if (nn .ne. n) call errquit('gasvdsol: b does not conform',nn,
645
1083
      nsing = min(m,n)
646
1084
      if (.not. ma_push_get(MT_DBL, nsing, 'gasvdsol', l_val, k_val))
647
1085
     $     call errquit('gasvdsol: val',nsing, MA_ERR)
648
 
      if (.not. ga_create(MT_DBL,m,nsing,'gasvd',0,0,g_u))
 
1086
      if (.not. ga_create(type,m,nsing,'gasvd',0,0,g_u))
649
1087
     $     call errquit('gasvdsol: u',m*nsing, GA_ERR)
650
 
      if (.not. ga_create(MT_DBL,nsing,n,'gasvd',0,0,g_vt))
 
1088
      if (.not. ga_create(type,nsing,n,'gasvd',0,0,g_vt))
651
1089
     $     call errquit('gasvdsol: u',nsing*n, GA_ERR)
652
 
      if (.not. ga_create(MT_DBL,nsing,nvec,'gasvd',0,0,g_tmp))
 
1090
      if (.not. ga_create(type,nsing,nvec,'gasvd',0,0,g_tmp))
653
1091
     $     call errquit('gasvdsol: tmp',nsing*nvec, GA_ERR)
654
1092
      call ga_zero(g_tmp)
655
1093
c
656
1094
      call ga_svd_seq(g_a, g_u,g_vt,dbl_mb(k_val))
657
1095
c
658
1096
      do i = 0, nsing-1
659
 
         if (dbl_mb(k_val+i) .lt. tol) then
660
 
            if (ga_nodeid() .eq. 0 .and. oprint) then
661
 
              write(6,*) ' neglecting ', i+1, dbl_mb(k_val+i)
662
 
            endif
663
 
            dbl_mb(k_val+i) = 0.0d0
664
 
         else
665
 
            dbl_mb(k_val+i) = 1.0d0/dbl_mb(k_val+i)
666
 
         end if
 
1097
        if (dbl_mb(k_val+i) .lt. tol) then
 
1098
          if (ga_nodeid() .eq. 0 .and. oprint) then
 
1099
            write(6,*) ' neglecting ', i+1, dbl_mb(k_val+i)
 
1100
          endif
 
1101
          dbl_mb(k_val+i) = 0.0d0
 
1102
        else
 
1103
          dbl_mb(k_val+i) = 1.0d0/dbl_mb(k_val+i)
 
1104
        end if
667
1105
      end do
668
1106
c
669
 
      call ga_dgemm('t','n',nsing,nvec,m,1.0d0,g_u,g_b,0.0d0,g_tmp)
 
1107
      call ga_dgemm('c','n',nsing,nvec,m,one,g_u,g_b,zero,g_tmp)
670
1108
      call ga_scale_lh(g_tmp,dbl_mb(k_val))
671
1109
      call ga_zero(g_x)
672
 
      call ga_dgemm('t','n',n,nvec,nsing,1.0d0,g_vt,g_tmp,0.0d0,g_x)
 
1110
      call ga_dgemm('c','n',n,nvec,nsing,one,g_vt,g_tmp,zero,g_x)
673
1111
c
674
1112
      if (.not. ga_destroy(g_tmp)) call errquit('gasvdsol: des',1,
675
1113
     &       GA_ERR)
681
1119
     &       GA_ERR)
682
1120
c
683
1121
      end
 
1122
C>
 
1123
C> \brief Perform SVD on rectangular matrix
 
1124
C>
684
1125
      subroutine ga_svd_seq(g_a, g_u, g_vt, values)
685
1126
      implicit none
686
1127
#include "errquit.fh"
704
1145
c
705
1146
      integer n, m, type, l_a, k_a, l_u, k_u, l_vt, k_vt, 
706
1147
     $     l_work, k_work, lwork, info, nsing
 
1148
      integer l_rwork, k_rwork
707
1149
c     
708
1150
      call ga_inquire(g_a, type, m, n)
709
1151
      nsing = min(m,n)
710
1152
      if (ga_nodeid() .eq. 0) then
711
1153
         lwork = 10*max(m,n)
712
 
         if (.not. ma_push_get(MT_DBL, m*n, 'gasvd', l_a, k_a))
 
1154
         if (.not. ma_push_get(type, m*n, 'gasvd', l_a, k_a))
713
1155
     $        call errquit('gasvd: a',m*n, MA_ERR)
714
 
         if (.not. ma_push_get(MT_DBL, m*nsing, 'gasvd', l_u, k_u))
 
1156
         if (.not. ma_push_get(type, m*nsing, 'gasvd', l_u, k_u))
715
1157
     $        call errquit('gasvd: u',m*nsing, MA_ERR)
716
 
         if (.not. ma_push_get(MT_DBL, nsing*n, 'gasvd', l_vt, k_vt))
 
1158
         if (.not. ma_push_get(type, nsing*n, 'gasvd', l_vt, k_vt))
717
1159
     $        call errquit('gasvd: vt',nsing*n, MA_ERR)
718
 
         if (.not. ma_push_get(MT_DBL, lwork, 'gasvd', l_work, k_work))
 
1160
         if (.not. ma_push_get(type, lwork, 'gasvd', l_work, k_work))
719
1161
     $        call errquit('gasvd: work',lwork, MA_ERR)
720
 
c
721
 
         call ga_get(g_a, 1, m, 1, n, dbl_mb(k_a), m)
722
 
c
723
 
         call dgesvd('s','s',m,n,dbl_mb(k_a),m,values,
724
 
     $        dbl_mb(k_u),m,dbl_mb(k_vt),nsing,
725
 
     $        dbl_mb(k_work),lwork,info)
726
 
         if (info .ne. 0) call errquit('gasvd: failed', info, MEM_ERR)
727
 
c
728
 
         call ga_put(g_u,  1, n,     1, nsing, dbl_mb(k_u),  n)
729
 
         call ga_put(g_vt, 1, nsing, 1, m,     dbl_mb(k_vt), n)
 
1162
         if (type.eq.MT_DCPL) then
 
1163
            if (.not. ma_push_get(MT_DBL, lwork, 'gasvd',
 
1164
     $                            l_rwork, k_rwork))
 
1165
     $           call errquit('gasvd: work',lwork, MA_ERR)
 
1166
         endif
 
1167
c
 
1168
         if (type.eq.MT_DBL) then
 
1169
c
 
1170
           call ga_get(g_a, 1, m, 1, n, dbl_mb(k_a), m)
 
1171
c
 
1172
           call dgesvd('s','s',m,n,dbl_mb(k_a),m,values,
 
1173
     $          dbl_mb(k_u),m,dbl_mb(k_vt),nsing,
 
1174
     $          dbl_mb(k_work),lwork,info)
 
1175
           if (info .ne. 0) call errquit('gasvd:d: failed',info,MEM_ERR)
 
1176
c
 
1177
           call ga_put(g_u,  1, n,     1, nsing, dbl_mb(k_u),  n)
 
1178
           call ga_put(g_vt, 1, nsing, 1, m,     dbl_mb(k_vt), n)
 
1179
c
 
1180
         else if (type.eq.MT_DCPL) then
 
1181
c
 
1182
           call ga_get(g_a, 1, m, 1, n, dcpl_mb(k_a), m)
 
1183
c
 
1184
           call zgesvd('s','s',m,n,dcpl_mb(k_a),m,values,
 
1185
     $          dcpl_mb(k_u),m,dcpl_mb(k_vt),nsing,
 
1186
     $          dcpl_mb(k_work),lwork,dbl_mb(k_rwork),info)
 
1187
           if (info .ne. 0) call errquit('gasvd:z: failed',info,MEM_ERR)
 
1188
c
 
1189
           call ga_put(g_u,  1, n,     1, nsing, dcpl_mb(k_u),  n)
 
1190
           call ga_put(g_vt, 1, nsing, 1, m,     dcpl_mb(k_vt), n)
 
1191
c
 
1192
         else
 
1193
c
 
1194
           call errquit('gasvd: illegal data type',type,UERR)
 
1195
c
 
1196
         endif
730
1197
c
731
1198
         if (.not. ma_chop_stack(l_a)) call errquit('gasvd ma',0,
732
1199
     &       MA_ERR)
736
1203
      call ga_sync()
737
1204
c     
738
1205
      end
 
1206
C> @}