← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwdft/rt_tddft/matutils/zmat_inv.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/nwdft/rt_tddft/matutils/zmat_inv.F
 
1
C====================================================================
 
2
C
 
3
C     Wrapper routine for inverting generic complex matrix g_za.
 
4
C
 
5
      subroutine zmat_inv (g_za, g_zainv)
 
6
      implicit none
 
7
 
 
8
#include "global.fh"
 
9
#include "errquit.fh"
 
10
#include "mafdecls.fh"
 
11
#include "stdio.fh"
 
12
#include "matutils.fh"
 
13
 
 
14
 
 
15
C     == Inputs ==
 
16
      integer, intent(in) :: g_za              !complex matrix to invert
 
17
 
 
18
 
 
19
C     == Outputs ==
 
20
      integer, intent(in) :: g_zainv           !inverse of matrix
 
21
 
 
22
 
 
23
C     == Parameter ==
 
24
      character(*), parameter :: pname = "zmat_inv: "
 
25
 
 
26
 
 
27
C     == Variables ==
 
28
      integer dim1, dim2, dtype
 
29
      integer n
 
30
 
 
31
 
 
32
C      
 
33
C     Get dims of GAs and check that they are correct types
 
34
C
 
35
C     Check the matrix (input 1).
 
36
C
 
37
      call ga_check_handle (g_za,
 
38
     $     "first argument of "//pname//"() is not a valid GA")
 
39
      call ga_inquire (g_za, dtype, dim1, dim2)
 
40
      if (dtype .ne. mt_dcpl) call errquit (pname//
 
41
     $     "expecting complex-valued GA as first argument", 0, 0)
 
42
      if (dim1 .ne. dim2)
 
43
     $     call errquit (pname//"dim1 must equal dim2", 0, 0)
 
44
 
 
45
C      
 
46
C     The size of all matricies must be n x n.
 
47
C
 
48
      n = dim1
 
49
 
 
50
      
 
51
C
 
52
C     Check the inverse (output) matrix.
 
53
C
 
54
      call ga_check_handle (g_zainv,
 
55
     $     "second argument of "//pname//"() is not a valid GA")
 
56
      call ga_inquire (g_zainv, dtype, dim1, dim2)
 
57
      if (dtype .ne. mt_dcpl) call errquit (pname//
 
58
     $     "expecting complex-valued GA as second argument", 0, 0)
 
59
      if (dim1 .ne. dim2)
 
60
     $     call errquit (pname//"dim1 must equal dim2", 0, 0)
 
61
      if (dim1.ne.n)
 
62
     $     call errquit (pname//"size of ainv must match size of a")
 
63
 
 
64
 
 
65
CXXX  [KAL]: only serial inverse routine works for now
 
66
      call zmat_inv_serial (n, g_za, g_zainv)
 
67
      
 
68
      end subroutine
 
69
 
 
70
 
 
71
C====================================================================
 
72
      subroutine zmat_inv_serial (n, g_za, g_zainv)
 
73
      implicit none
 
74
 
 
75
#include "global.fh"
 
76
#include "errquit.fh"
 
77
#include "mafdecls.fh"
 
78
#include "stdio.fh"
 
79
#include "matutils.fh"
 
80
 
 
81
 
 
82
C     == Inputs ==
 
83
      integer, intent(in) :: n                 !size of mats
 
84
      integer, intent(in) :: g_za              !input matrix
 
85
 
 
86
 
 
87
C     == Outputs ==      
 
88
      integer, intent(in) :: g_zainv           !inverse of matrix
 
89
 
 
90
 
 
91
C     == Parameter ==
 
92
      character(*), parameter :: pname = "zmat_inv_serial: "
 
93
 
 
94
 
 
95
C     == Variables ==
 
96
      logical ok
 
97
      integer i, j
 
98
      integer me
 
99
      integer info
 
100
      integer g_zident
 
101
      integer izident, lzident
 
102
      integer lza, iza
 
103
      integer lzainv, izainv
 
104
      integer lpivot, ipivot
 
105
 
 
106
      
 
107
      call ga_sync ()
 
108
      
 
109
C     
 
110
C     Do all work on all processors (wasteful but works).
 
111
C     
 
112
      me = ga_nodeid ()
 
113
 
 
114
 
 
115
C     
 
116
C     Allocation.
 
117
C      
 
118
      ok=.true.
 
119
      ok=ok.and.ma_push_get(mt_dcpl, n*n, "complex mat A", lza,iza)
 
120
      ok=ok.and.ma_push_get(mt_dcpl, n*n, "complex mat Ainv",
 
121
     $     lzainv, izainv)
 
122
      ok=ok.and.ma_push_get(mt_dcpl, n*n, "complex ident mat",
 
123
     $     lzident, izident)
 
124
      ok=ok.and.ma_push_get(mt_int, n, "pivot indices", lpivot, ipivot)
 
125
      
 
126
      ok = ok .and. ga_duplicate (g_za, g_zident, "ident")
 
127
      if (.not.ok) call errquit (pname//"alloc failed", 0, 0)
 
128
      
 
129
      
 
130
C
 
131
C     Unpack buffer to GA
 
132
C
 
133
      call pack_ga2buffer_dcpl (g_za, dcpl_mb(iza))
 
134
 
 
135
      
 
136
C
 
137
C     LU factorization
 
138
C
 
139
      call zgetrf (n, n, dcpl_mb(iza), n, int_mb(ipivot), info)
 
140
 
 
141
      if (info .ne. 0)
 
142
     $     call errquit (pname//"LU factorization failed", 0, 0)
 
143
 
 
144
 
 
145
C
 
146
C     Make complex identity matrix buffer.
 
147
C
 
148
      call ga_zero (g_zident)
 
149
      call mat_set_ident (g_zident)
 
150
      call pack_ga2buffer_dcpl (g_zident, dcpl_mb(izident))
 
151
 
 
152
      
 
153
C
 
154
C     Compute inverse by solving AX = I (I is the nxn ident mat)
 
155
C
 
156
C      SUBROUTINE ZGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
 
157
 
 
158
      call zgetrs ("N", n, n, dcpl_mb(iza), n, int_mb(ipivot),
 
159
     $     dcpl_mb(izident), n, info)
 
160
      
 
161
      if (info .ne. 0)
 
162
     $     call errquit (pname//"Inversion failed", 0, 0)
 
163
 
 
164
 
 
165
C
 
166
C     Pack solution (stored in zident buffer) into output GA.
 
167
C
 
168
      call pack_buffer2ga_dcpl (dcpl_mb(izident), g_zainv)
 
169
 
 
170
 
 
171
C
 
172
C     Clean up.
 
173
C      
 
174
      ok = .true.
 
175
      ok = ok .and. ma_chop_stack (lza)
 
176
      ok = ok .and. ga_destroy (g_zident)
 
177
      if (.not.ok) call errquit (pname//"free failed", 0, 0)
 
178
      
 
179
      end subroutine
 
180
 
 
181
 
 
182
C====================================================================
 
183
      subroutine zmat_inv_check (g_za, g_zainv)
 
184
      implicit none
 
185
 
 
186
#include "global.fh"
 
187
#include "errquit.fh"
 
188
#include "mafdecls.fh"
 
189
#include "stdio.fh"
 
190
#include "matutils.fh"
 
191
 
 
192
C     == Inputs ==
 
193
      integer, intent(in) :: g_za              !complex matrix to invert
 
194
      integer, intent(in) :: g_zainv           !inverse of matrix
 
195
 
 
196
 
 
197
C     == Parameter ==
 
198
      character(*), parameter     :: pname = "zmat_inv_check: "
 
199
      double complex, parameter   :: z1 = (1d0, 0d0)
 
200
      double complex, parameter   :: z0 = (0d0, 0d0)
 
201
      double precision, parameter :: thresh = 1d-8
 
202
 
 
203
 
 
204
C     == Variables ==
 
205
      integer dim1, dim2, dtype
 
206
      integer n
 
207
      integer g_zprod
 
208
 
 
209
 
 
210
C      
 
211
C     Get dims of GAs and check that they are correct types
 
212
C
 
213
C     Check the matrix (input 1).
 
214
C
 
215
      call ga_check_handle (g_za,
 
216
     $     "first argument of "//pname//"() is not a valid GA")
 
217
      call ga_inquire (g_za, dtype, dim1, dim2)
 
218
      if (dtype .ne. mt_dcpl) call errquit (pname//
 
219
     $     "expecting complex-valued GA as first argument", 0, 0)
 
220
      if (dim1 .ne. dim2)
 
221
     $     call errquit (pname//"dim1 must equal dim2", 0, 0)
 
222
 
 
223
C      
 
224
C     The size of all matricies must be n x n.
 
225
C
 
226
      n = dim1
 
227
 
 
228
      
 
229
C
 
230
C     Check the inverse (output) matrix.
 
231
C
 
232
      call ga_check_handle (g_zainv,
 
233
     $     "second argument of "//pname//"() is not a valid GA")
 
234
      call ga_inquire (g_zainv, dtype, dim1, dim2)
 
235
      if (dtype .ne. mt_dcpl) call errquit (pname//
 
236
     $     "expecting complex-valued GA as second argument", 0, 0)
 
237
      if (dim1 .ne. dim2)
 
238
     $     call errquit (pname//"dim1 must equal dim2", 0, 0)
 
239
      if (dim1.ne.n)
 
240
     $     call errquit (pname//"size of ainv must match size of a")
 
241
 
 
242
 
 
243
C
 
244
C     Allocation
 
245
C
 
246
      if (.not. ga_duplicate (g_za, g_zprod, "prod"))
 
247
     $     call errquit (pname//"failed to alloc prod", 0, GA_ERR)
 
248
 
 
249
 
 
250
 
 
251
C
 
252
C     Check that A^-1 A = I
 
253
C      
 
254
      call ga_zgemm ("N", "N", n, n, n, z1, g_zainv, g_za, z0, g_zprod)
 
255
      if (.not. mat_is_ident (g_zprod, thresh))
 
256
     $     call errquit (pname//"inversion check failed", 0, 0)
 
257
      
 
258
 
 
259
C
 
260
C     Clean up
 
261
C     
 
262
      if (.not. ga_destroy (g_zprod))
 
263
     $     call errquit (pname//"destroy failed", 0, 0)
 
264
 
 
265
      end subroutine
 
266
 
 
267
 
 
268
C====================================================================
 
269
      subroutine zmat_inv_example_driver ()
 
270
      implicit none
 
271
 
 
272
#include "global.fh"
 
273
#include "errquit.fh"
 
274
#include "mafdecls.fh"
 
275
#include "stdio.fh"
 
276
 
 
277
 
 
278
C     == Parameters ==
 
279
      character(*), parameter :: pname = "zmat_inv_example_driver: "
 
280
 
 
281
 
 
282
C     == Variables ==
 
283
      integer me
 
284
      integer g_zmat, g_zmatinv
 
285
      double complex val
 
286
 
 
287
 
 
288
      me = ga_nodeid ()
 
289
 
 
290
 
 
291
C
 
292
C     EXAMPLE ONE
 
293
C
 
294
C     octave:1> A=[[1,2,0];[0,3,0];[2,-4,2]]
 
295
C     A =
 
296
C
 
297
C     1   2   0
 
298
C     0   3   0
 
299
C     2  -4   2
 
300
C
 
301
C     octave:2> inverse(A)
 
302
C     ans =
 
303
C
 
304
C     1.00000  -0.66667   0.00000
 
305
C     0.00000   0.33333   0.00000
 
306
C     -1.00000   1.33333   0.50000
 
307
C
 
308
 
 
309
      if (me.eq.0) then
 
310
         write (luout, *) ""
 
311
         write (luout, *) ""
 
312
         call util_print_centered (luout, "Inversion example one",
 
313
     $        20, .true.)
 
314
      endif
 
315
      
 
316
 
 
317
C
 
318
C     Allocation.
 
319
C
 
320
      if (.not. ga_create(mt_dcpl,3,3,"matrix" , 0, 0, g_zmat))
 
321
     $     call errquit ("failed to create mat", 0, 0)
 
322
 
 
323
      if (.not. ga_create(mt_dcpl,3,3,"matrix inv" , 0, 0, g_zmatinv))
 
324
     $     call errquit ("failed to create mat inv", 0, 0)
 
325
      
 
326
      call ga_zero (g_zmat)
 
327
 
 
328
      val = dcmplx (1d0, 0d0)
 
329
      call ga_put (g_zmat, 1, 1, 1, 1, val, 1)
 
330
 
 
331
      val = dcmplx (2d0, 0d0)
 
332
      call ga_put (g_zmat, 1, 1, 2, 2, val, 1)
 
333
 
 
334
      val = dcmplx (3d0, 0d0)
 
335
      call ga_put (g_zmat, 2, 2, 2, 2, val, 1)
 
336
 
 
337
      val = dcmplx (2d0, 0d0)
 
338
      call ga_put (g_zmat, 3, 3, 1, 1, val, 1)
 
339
 
 
340
      val = dcmplx (-4d0, 0d0)
 
341
      call ga_put (g_zmat, 3, 3, 2, 2, val, 1)
 
342
 
 
343
      val = dcmplx (2d0, 0d0)
 
344
      call ga_put (g_zmat, 3, 3, 3, 3, val, 1)
 
345
 
 
346
 
 
347
C
 
348
C     Inversion
 
349
C
 
350
      call zmat_inv (g_zmat, g_zmatinv)
 
351
      call zmat_inv_check (g_zmat, g_zmatinv)
 
352
      call ga_print (g_zmat)
 
353
      call ga_print (g_zmatinv)
 
354
 
 
355
 
 
356
C
 
357
C     EXAMPLE TWO
 
358
C
 
359
C     octave:1> A=[[1.2+i,1,i];[0,3,-0.98*i];[2+3.2*i,-4,2.5+8*i]]
 
360
C     A =
 
361
C
 
362
C     1.20000 + 1.00000i   1.00000 + 0.00000i   0.00000 + 1.00000i
 
363
C     0.00000 + 0.00000i   3.00000 + 0.00000i  -0.00000 - 0.98000i
 
364
C     2.00000 + 3.20000i  -4.00000 + 0.00000i   2.50000 + 8.00000i
 
365
C     
 
366
C     octave:2> inverse(A)
 
367
C     ans =
 
368
C     
 
369
C     0.867525 - 0.256532i  -0.512594 + 0.069857i  -0.167565 - 0.011740i
 
370
C     -0.073269 - 0.137812i   0.419917 + 0.105576i   0.046621 + 0.044729i
 
371
C     -0.421875 + 0.224292i   0.323190 - 0.265053i   0.136924 - 0.142717i
 
372
C
 
373
 
 
374
      if (me.eq.0) then
 
375
         write (luout, *) ""
 
376
         write (luout, *) ""
 
377
         call util_print_centered (luout, "Inversion example two",
 
378
     $        20, .true.)
 
379
      endif
 
380
 
 
381
      call ga_zero (g_zmat)
 
382
 
 
383
      val = dcmplx (1.2d0, 1d0)
 
384
      call ga_put (g_zmat, 1, 1, 1, 1, val, 1)
 
385
 
 
386
      val = dcmplx (1.0d0, 0d0)
 
387
      call ga_put (g_zmat, 1, 1, 2, 2, val, 1)
 
388
 
 
389
      val = dcmplx (0d0, 1d0)
 
390
      call ga_put (g_zmat, 1, 1, 3, 3, val, 1)
 
391
 
 
392
      val = dcmplx (0d0, 0d0)
 
393
      call ga_put (g_zmat, 2, 2, 1, 1, val, 1)
 
394
 
 
395
      val = dcmplx (3d0, 0d0)
 
396
      call ga_put (g_zmat, 2, 2, 2, 2, val, 1)
 
397
 
 
398
      val = dcmplx (0d0, -0.98d0)
 
399
      call ga_put (g_zmat, 2, 2, 3, 3, val, 1)
 
400
 
 
401
      val = dcmplx (2d0, 3.2d0)
 
402
      call ga_put (g_zmat, 3, 3, 1, 1, val, 1)
 
403
 
 
404
      val = dcmplx (-4d0, 0d0)
 
405
      call ga_put (g_zmat, 3, 3, 2, 2, val, 1)
 
406
 
 
407
      val = dcmplx (2.5d0, 8d0)
 
408
      call ga_put (g_zmat, 3, 3, 3, 3, val, 1)
 
409
 
 
410
 
 
411
C
 
412
C     Inversion
 
413
C
 
414
      call zmat_inv (g_zmat, g_zmatinv)
 
415
      call zmat_inv_check (g_zmat, g_zmatinv)
 
416
      call ga_print (g_zmat)
 
417
      call ga_print (g_zmatinv)
 
418
 
 
419
 
 
420
C
 
421
C     Clean up.
 
422
C     
 
423
      if (.not. ga_destroy (g_zmat))
 
424
     $     call errquit (pname//"failed to destroy g_zmat",0,0)
 
425
      if (.not. ga_destroy (g_zmatinv))
 
426
     $     call errquit (pname//"failed to destroy g_zmatinv",0,0)
 
427
 
 
428
      end subroutine
 
429
      
 
430