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

« back to all changes in this revision

Viewing changes to src/tools/ga-5-1/armci/examples/benchmarks/cg/armci_sharedmemory/cg.c

  • 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
* 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:
1
 
#if HAVE_CONFIG_H
2
 
#   include "config.h"
3
 
#endif
4
 
 
5
 
/** @file
6
 
 * $Id: cg.c,v 1.1.2.2 2007-07-05 15:29:56 vinod Exp $
7
 
 */
8
 
 
9
 
#if HAVE_STDIO_H
10
 
#   include <stdio.h>
11
 
#endif
12
 
#if HAVE_MATH_H
13
 
#   include <math.h>
14
 
#endif
15
 
#if HAVE_STDLIB_H
16
 
#   include <stdlib.h>
17
 
#endif
18
 
#if HAVE_SYS_TYPES_H
19
 
#   include <sys/types.h>
20
 
#endif
21
 
#if HAVE_SYS_STAT_H
22
 
#   include <sys/stat.h>
23
 
#endif
24
 
#if HAVE_FCNTL_H
25
 
#   include <fcntl.h>
26
 
#endif
27
 
#if HAVE_STRING_H
28
 
#   include <string.h>
29
 
#endif
30
 
 
31
 
#include "armci.h"
32
 
#include "message.h"
33
 
 
34
 
#define PRINT_VEC_
35
 
 
36
 
int na,nz;
37
 
double *bvec,*dvec,*svec,*dmvec,*m_dvec,*amat,*xvec,*axvec,*rvec,*qvec;
38
 
int *ridx,*cidx;
39
 
int me, nproc;
40
 
int myfirstrow=0,mylastrow=0;
41
 
double epsilon=1e-4;
42
 
double time_get=0;
43
 
static int niter;
44
 
void read_and_create(int,char **);
45
 
void computeminverser(double *,double *, double *);
46
 
void computeminverse(double *,double *, int *,int *);
47
 
void finalize_arrays(int);
48
 
extern void acg_matvecmul(double *,double *,double *,int *,int *);
49
 
extern void acg_addvec(double *,double *,double *,double *, double *);
50
 
extern void acg_2addvec(double *,double *, double *,double *, double *,
51
 
        double *, double *,double *, double *, double *, int *, int *);
52
 
extern double acg_ddot(double *,double *);
53
 
 
54
 
 
55
 
void conjugate_gradient(int nit,int dopreconditioning)
56
 
{
57
 
    int i;
58
 
    double d_one=1.0,d_negone=-1.0;
59
 
    double delta0=0.0,deltaold=0.0,deltanew=0.0;
60
 
    double alpha=0.0,negalpha,beta,dtransposeq;
61
 
 
62
 
    acg_matvecmul(amat,xvec,axvec,ridx,cidx);             /* compute Ax */
63
 
#ifdef PRINT_VEC 
64
 
    acg_printvec("axvec",axvec);
65
 
    acg_printvec("bvec",bvec);
66
 
#endif
67
 
 
68
 
    /* TODO original call had too many arguments */
69
 
    /* acg_addvec(&d_one,bvec,&d_negone,axvec,rvec,ridx,cidx); */ /* r=b-Ax */
70
 
    acg_addvec(&d_one,bvec,&d_negone,axvec,rvec); /* r=b-Ax */
71
 
 
72
 
#ifdef PRINT_VEC 
73
 
    acg_printvec("rvec",rvec);
74
 
    if(me==0)for(i=0;i<nz;i++)printf("\n%d:col[%d]=%d",me,i,cidx[i]);
75
 
    fflush(stdout);
76
 
#endif
77
 
 
78
 
    if(dopreconditioning){
79
 
        computeminverse(dmvec,amat,ridx,cidx);
80
 
        /*acg_printvec("dmvec",dmvec);*/
81
 
        computeminverser(dmvec,rvec,dvec);
82
 
        /*acg_printvec("dvec",dvec);*/
83
 
        if (me == 0)
84
 
            printf("\nDoing preconditioning!\n");
85
 
    }
86
 
    else{
87
 
        if(me==0)memcpy(dvec,rvec,na*sizeof(double));
88
 
    }
89
 
    if(dopreconditioning)
90
 
        deltanew = acg_ddot(rvec,dvec); /* deltanew = r.r_tranpose */
91
 
    else
92
 
        deltanew = acg_ddot(rvec,rvec); /* deltanew = r.r_tranpose */
93
 
 
94
 
    delta0 = deltanew;                  /* delta0 = deltanew */
95
 
 
96
 
    if(me==0)printf("\n\tdelta0 is %f\n",delta0);
97
 
 
98
 
    for(i=0;i<nit && deltanew>(epsilon*epsilon*delta0);i++){
99
 
        acg_matvecmul(amat,dvec,qvec,ridx,cidx); /* q = ad */
100
 
 
101
 
        dtransposeq=acg_ddot(dvec,qvec);         /* compute d_transpose.q */
102
 
 
103
 
        alpha = deltanew/dtransposeq;            /* deltanew/(d_transpose.q) */
104
 
#if 0
105
 
        if(i>0 && i%50==0){
106
 
            /* compute Ax*/
107
 
            acg_matvecmul(amat,xvec,axvec,ridx,cidx);
108
 
            /* x = x+ alpha.d*/ /* r=b-Ax*/
109
 
            acg_2addvec(&d_one,xvec,&alpha,dvec,xvec,&d_one,bvec,
110
 
                    &d_negone,axvec,rvec,ridx,cidx);
111
 
        }
112
 
        else
113
 
#endif
114
 
        {
115
 
            negalpha = 0.0-alpha;                         
116
 
            /* x = x+ alpha.d*/ /* r=r-alpha.q*/
117
 
            acg_2addvec(&d_one,xvec,&alpha,dvec,xvec,&d_one,rvec,
118
 
                    &negalpha,qvec,rvec,ridx,cidx);
119
 
        }
120
 
 
121
 
        if(dopreconditioning)
122
 
            computeminverser(dmvec,rvec,svec);
123
 
 
124
 
        deltaold = deltanew;                /* deltaold = deltanew*/
125
 
 
126
 
        if(dopreconditioning)
127
 
            deltanew = acg_ddot(svec,rvec); /* deltanew = r_transpose.r*/
128
 
        else
129
 
            deltanew = acg_ddot(rvec,rvec); /* deltanew = r_transpose.r*/
130
 
 
131
 
        beta = deltanew/deltaold;           /* beta = deltanew/deltaold*/
132
 
 
133
 
        if(dopreconditioning) {
134
 
            /* too many aguments in function call */
135
 
            /* acg_addvec(&d_one,svec,&beta,dvec,dvec,ridx,cidx); */
136
 
            /* d = s + beta.d */
137
 
            acg_addvec(&d_one,svec,&beta,dvec,dvec);
138
 
        } else {
139
 
            /* too many aguments in function call */
140
 
            /* acg_addvec(&d_one,rvec,&beta,dvec,dvec,ridx,cidx); */
141
 
            /* d = r + beta.d */
142
 
            acg_addvec(&d_one,rvec,&beta,dvec,dvec);
143
 
        }
144
 
 
145
 
#ifdef PRINT_VEC 
146
 
        acg_printvec("xvec",xvec);
147
 
#endif
148
 
        /* acg_printvec("xvec",xvec); */
149
 
 
150
 
    }
151
 
    if(me==0)printf("\n\tIteration:%d\tBeta:%0.4f\tDelta:%f\n",
152
 
            i,beta,deltanew);
153
 
    niter = i;
154
 
}
155
 
 
156
 
 
157
 
void initialize_arrays(int dpc)
158
 
{
159
 
    int i;
160
 
    for(i=0;i<na;i++){
161
 
        xvec[i]=0.0;
162
 
        dvec[i]=axvec[i]=rvec[i]=qvec[i]=0;
163
 
        if(dpc){dmvec[i]=svec[i]=0;}
164
 
    }
165
 
}
166
 
 
167
 
 
168
 
void finalize_arrays(int dpc)
169
 
{
170
 
    if(me==0){
171
 
        ARMCI_Free(bvec);
172
 
        ARMCI_Free(dvec);
173
 
        ARMCI_Free(amat);
174
 
        ARMCI_Free(xvec);
175
 
        ARMCI_Free(axvec);
176
 
        ARMCI_Free(rvec);
177
 
        ARMCI_Free(qvec);
178
 
        ARMCI_Free(ridx);
179
 
        ARMCI_Free(cidx);
180
 
        if(dpc){
181
 
            ARMCI_Free(svec);
182
 
            ARMCI_Free(dmvec);
183
 
        }
184
 
    }
185
 
}     
186
 
 
187
 
 
188
 
FILE *fd;
189
 
 
190
 
 
191
 
int main(int argc, char **argv)
192
 
{
193
 
    int dopreconditioning=1;
194
 
    double time0,time1;
195
 
 
196
 
    armci_msg_init(&argc,&argv);
197
 
    nproc = armci_msg_nproc();
198
 
    me = armci_msg_me();
199
 
    ARMCI_Init(); /* initialize ARMCI */
200
 
 
201
 
    if(me==0)printf("\n                          CONJUGATE GRADIENT EXAMPLE\n");
202
 
    if(argc<3){
203
 
        if(me==0){
204
 
            printf(" CORRECT USAGE IS:");
205
 
            printf("\n\n <launch commands> cg.x na nz file");
206
 
            printf("\n\n where:");
207
 
            printf("\n\tna is array dimention (only square arrays supported)");
208
 
            printf("\n\tnz is number of non-zeros");
209
 
            printf("\n\tfile is either the input file or the word random");
210
 
            printf("\n\t  use the word random if you to use random input");
211
 
            printf("\n\t  input should be in row compressed format");
212
 
            printf("\n\t  file should have matrix a followed by row, col & b (Ax=b)");
213
 
            printf("\n\t  if file also has na and nz, pass them as 0's and the");
214
 
            printf("\n\t  program will read them from the file");
215
 
            printf("\n\nexample usages are:");
216
 
            printf("\n\tmpirun -np 4 ./ga_cg.x 5000 80000 /home/me/myinput.dat");
217
 
            printf("\n\tor");
218
 
            printf("\n\tmpirun -np 4 ./ga_cg.x 5000 80000 random\n\n");
219
 
            fflush(stdout);
220
 
        }
221
 
        ARMCI_Finalize();
222
 
        armci_msg_finalize();
223
 
        return 0;
224
 
    }
225
 
 
226
 
    read_and_create(argc,argv);
227
 
 
228
 
    if(me==0)printf("\nWarmup and initialization run");
229
 
#if 0
230
 
    initialize_arrays(dopreconditioning);
231
 
    conjugate_gradient(1,dopreconditioning);
232
 
    time_get =0.0;
233
 
#endif
234
 
    if(me==0)printf("\n\nStarting Conjugate Gradient ....");
235
 
    initialize_arrays(dopreconditioning);
236
 
    time0=armci_timer();
237
 
    conjugate_gradient(30000/*2*/,dopreconditioning);
238
 
    time1=armci_timer();
239
 
 
240
 
    acg_matvecmul(amat,xvec,axvec,ridx,cidx);
241
 
    if(me==0)printf("\n%d:in %d iterations time to solution=%f-%f ax and b in cg_output.out\n",me,niter,(time1-time0),time_get);
242
 
    acg_matvecmul(amat,xvec,axvec,ridx,cidx);
243
 
    if(me==0){
244
 
        int i;
245
 
        fd = fopen("cg_output.out", "w");
246
 
        for(i=0;i<=na;i++)
247
 
            fprintf(fd,"\n%d:%s[%d]=%f %s[%d]=%f",me,"bvec",i,bvec[i],"axvec",i,axvec[i]);
248
 
        fflush(stdout);
249
 
        fclose(fd);
250
 
    }
251
 
 
252
 
    finalize_arrays(dopreconditioning);
253
 
    armci_msg_barrier();
254
 
 
255
 
    if(me==0)printf("Terminating ..\n");
256
 
    ARMCI_Finalize();
257
 
    armci_msg_finalize();
258
 
    return 0;
259
 
}