← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/tools/ga-4-3/global/examples/petsc.F

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

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

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

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tools/ga-4-3/global/examples/petsc.F
 
1
      program petsc
 
2
      implicit none
 
3
#include "mafdecls.fh"
 
4
#include "global.fh"
 
5
 
 
6
      integer heap, stack
 
7
      data heap,stack /2*40000/
 
8
c
 
9
c***  Intitialize a message passing library
 
10
c
 
11
#ifdef MPI
 
12
      call mpi_init
 
13
#else
 
14
      call pbeginf
 
15
#endif
 
16
c
 
17
      call ga_initialize()                               ! initialize GAs 
 
18
      if (.not. ma_init(MT_DBL, heap, stack))            ! initialize MA
 
19
     $    call ga_error("ma init failed",heap+stack)
 
20
#ifdef GA_TRACE
 
21
      call trace_init(10000)                             ! initialize trace
 
22
#endif
 
23
      call iterate()                                     ! do the work
 
24
#ifdef GA_TRACE
 
25
      call trace_end(ga_nodeid())                        ! end trace
 
26
#endif
 
27
      call ga_terminate()                                ! terminate GAs 
 
28
c
 
29
#ifdef MPI
 
30
      call mpi_finalize()
 
31
#else
 
32
      call pend()
 
33
#endif
 
34
      end
 
35
 
 
36
      subroutine iterate()
 
37
      implicit none
 
38
 
 
39
#include "include/finclude/petsc.h"
 
40
#include "include/finclude/vec.h"
 
41
#include "include/finclude/mat.h"
 
42
#include "include/finclude/pc.h"
 
43
#include "include/finclude/ksp.h"
 
44
#include "include/finclude/sles.h"
 
45
#include "include/finclude/sys.h"
 
46
 
 
47
#include "mafdecls.fh"
 
48
#include "global.fh"
 
49
 
 
50
      double precision  norm
 
51
      integer     i, j, II, JJ, ierr, m, n
 
52
      parameter (m = 3)
 
53
      parameter (n = 3)
 
54
      integer     its, Istart, Iend, flg
 
55
      integer     me, nproc
 
56
      Scalar      v, one, neg_one
 
57
      Vec         x, b, u
 
58
      Mat         A 
 
59
      SLES        sles
 
60
      KSP         ksp
 
61
      PetscRandom rctx
 
62
 
 
63
      integer g_x
 
64
      integer ld
 
65
      double precision buf_v(1)
 
66
      PetscOffset idx
 
67
 
 
68
c$$$      PC          pc
 
69
c$$$      PCType      ptype 
 
70
c$$$      double precision tol
 
71
 
 
72
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
73
!                 Beginning of program
 
74
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
75
c***  check parallel environment
 
76
      me = ga_nodeid()
 
77
      nproc = ga_nnodes()
 
78
 
 
79
      one  = 1.0
 
80
      neg_one = -1.0
 
81
      ld = m*n
 
82
 
 
83
c***  Global Array
 
84
c***  create global arrays: g_x - approx. solution
 
85
      if (.not. ga_create(MT_DBL, m*n, 1, 'x', 1, 1, g_x))
 
86
     $     call ga_error(' ga_create failed ',0)
 
87
c
 
88
c***  initial guess for x -- zero
 
89
      call ga_zero(g_x)
 
90
c
 
91
c$$$      do i=1,m*n
 
92
c$$$         buf(i) = i
 
93
c$$$      enddo
 
94
c$$$      call ga_put(g_x,1,m*n,1,1,buf,ld)
 
95
c
 
96
c***  PETSC      
 
97
      call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 
98
 
 
99
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
100
!      Compute the matrix and right-hand-side vector that define
 
101
!      the linear system, Ax = b.
 
102
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
103
 
 
104
c***  Matrix A
 
105
      call MatCreate(PETSC_COMM_WORLD,m*n,m*n,A,ierr)
 
106
      call MatGetOwnershipRange(A,Istart,Iend,ierr)
 
107
 
 
108
      do II=Istart,Iend-1
 
109
        v = -1.0
 
110
        i = II/n
 
111
        j = II - i*n  
 
112
        if ( i.gt.0 ) then
 
113
          JJ = II - n
 
114
          call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr)
 
115
        endif
 
116
        if ( i.lt.m-1 ) then
 
117
          JJ = II + n
 
118
          call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr)
 
119
        endif
 
120
        if ( j.gt.0 ) then
 
121
          JJ = II - 1
 
122
          call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr)
 
123
        endif
 
124
        if ( j.lt.n-1 ) then
 
125
          JJ = II + 1
 
126
          call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr)
 
127
        endif
 
128
        v = 4.0
 
129
        call  MatSetValues(A,1,II,1,II,v,ADD_VALUES,ierr)
 
130
      enddo
 
131
 
 
132
      call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 
133
      call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
 
134
 
 
135
c***  Vector b
 
136
      call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr)
 
137
      call VecSetFromOptions(u,ierr)
 
138
      call VecDuplicate(u,b,ierr)
 
139
      call VecDuplicate(b,x,ierr)
 
140
 
 
141
c***  u is the exact solution
 
142
      call VecSet(one,u,ierr)
 
143
      if (me .eq. 0) print *, 'Exact solution:'
 
144
      call VecView(u, VIEWER_STDOUT_WORLD, ierr)
 
145
 
 
146
c***  b is the right hand side
 
147
      call MatMult(A,u,b,ierr)
 
148
 
 
149
c***  Manage to make connection of ga to petsc: g_x -> x
 
150
      call VecGetOwnershipRange(x,Istart,Iend,ierr)
 
151
      call VecGetArray(x,buf_v,idx,ierr)
 
152
 
 
153
      call ga_get(g_x,Istart+1,Iend,1,1,buf_v(idx+1),ld)
 
154
 
 
155
      call VecRestoreArray(x,buf_v,idx,ierr)
 
156
 
 
157
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
158
!         Create the linear solver and set various options
 
159
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
160
 
 
161
      call SLESCreate(PETSC_COMM_WORLD,sles,ierr)
 
162
 
 
163
      call SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN,
 
164
     &                      ierr)
 
165
 
 
166
c$$$      call SLESGetPC(sles,pc,ierr)
 
167
c$$$      ptype = PCJACOBI
 
168
c$$$      call PCSetType(pc,ptype,ierr)
 
169
c$$$      call SLESGetKSP(sles,ksp,ierr)
 
170
c$$$      tol = 1.e-7
 
171
c$$$      call KSPSetTolerances(ksp,tol,PETSC_DEFAULT_DOUBLE_PRECISION,
 
172
c$$$     &     PETSC_DEFAULT_DOUBLE_PRECISION,PETSC_DEFAULT_INTEGER,ierr)
 
173
 
 
174
      call SLESSetFromOptions(sles,ierr)
 
175
 
 
176
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
177
!                      Solve the linear system
 
178
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
179
 
 
180
      call SLESSolve(sles,b,x,its,ierr)
 
181
 
 
182
      call ga_sync()
 
183
c***  write the approx solution back to ga
 
184
      call VecGetArray(x,buf_v,idx,ierr)
 
185
 
 
186
      call ga_put(g_x,Istart+1,Iend,1,1,buf_v(idx+1),ld)
 
187
 
 
188
      call VecRestoreArray(x,buf_v,idx,ierr)
 
189
      
 
190
      if (me .eq. 0) print *, 'Approx solution:'
 
191
      call ga_print(g_x)
 
192
 
 
193
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
194
!                     Check solution and clean up
 
195
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
196
 
 
197
!  Check the error
 
198
 
 
199
      call VecAXPY(neg_one,u,x,ierr)
 
200
      call VecNorm(x,NORM_2,norm,ierr)
 
201
      if (me .eq. 0) then
 
202
        if (norm .gt. 1.e-12) then
 
203
           write(6,100) norm, its
 
204
        else
 
205
           write(6,110) its
 
206
        endif
 
207
      endif
 
208
  100 format('Norm of error ',e10.4,' iterations ',i5)
 
209
  110 format('Norm of error < 1.e-12, iterations ',i5)
 
210
 
 
211
c***  clean up
 
212
      call SLESDestroy(sles,ierr)
 
213
      call VecDestroy(u,ierr)
 
214
      call VecDestroy(x,ierr)
 
215
      call VecDestroy(b,ierr)
 
216
      call MatDestroy(A,ierr)
 
217
 
 
218
      call PetscFinalize(ierr)
 
219
 
 
220
      if(.not. ga_destroy(g_x)) call ga_error('invalid handle ?',0)
 
221
c
 
222
      call ga_sync()
 
223
      end
 
224