source: branches/MPI/ADOL-C/src/adolc_mpi.cpp @ 424

Last change on this file since 424 was 424, checked in by kulshres, 9 years ago

update with 'mpi' branch from 'gitclone'

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

File size: 16.7 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     adolc_mpi.cpp
4 Revision: $Id$
5 Contents: C allocation of arrays of doubles in several dimensions
6
7 Copyright (c) Andrea Walther, Benjamin Letschert
8
9 This file is part of ADOL-C. This software is provided as open source.
10 Any use, reproduction, or distribution of the software constitutes
11 recipient's acceptance of the terms of the accompanying license file.
12
13----------------------------------------------------------------------------*/
14
15#include <adolc/common.h>
16#include <adolc/adolc_mpi.h>
17#include "oplate.h"
18#include "taping_p.h"
19#include <adolc/adouble.h>
20#include <adolc/drivers/drivers.h>
21#include <adolc/tapedoc/tapedoc.h>
22#include <adolc/adalloc.h>
23#include <adolc/interfaces_mpi.h>
24#include <adolc/convolut.h>
25
26int mpi_initialized = 0;
27int all_root = 0;
28
29int trace_on( int id,
30              int size,
31              short tag,
32              int keepTaylors
33){
34    return trace_on(id+size*tag,keepTaylors);
35}
36
37BEGIN_C_DECLS
38
39int ADOLC_MPI_Init( int* a,
40                    char*** b
41){
42    mpi_initialized = 1;
43    all_root = 0;
44    return MPI_Init(a,b);
45}
46int ADOLC_MPI_Comm_size( ADOLC_MPI_Comm comm,
47                         int* size
48){
49    int ierr = MPI_Comm_size(comm,size);
50    return ierr;
51}
52int ADOLC_MPI_Comm_rank( ADOLC_MPI_Comm comm,
53                         int* rank
54){
55    return MPI_Comm_rank(comm, rank);
56}
57
58int ADOLC_MPI_Get_processor_name( char* a,
59                                  int* b
60){
61    return MPI_Get_processor_name(a,b);
62}
63
64int ADOLC_MPI_Barrier( ADOLC_MPI_Comm comm ){
65    put_op(barrier_op);
66    return MPI_Barrier(comm);
67}
68
69int ADOLC_MPI_Finalize( ){
70    return MPI_Finalize();
71}
72
73END_C_DECLS
74
75int ADOLC_MPI_Send( adouble *buf,
76                    int count,
77                    ADOLC_MPI_Datatype datatype,
78                    int dest,
79                    int tag,
80                    ADOLC_MPI_Comm comm
81){
82    int i,h=count;
83    int ierr =0;
84    double *trade;
85
86    trade = myalloc1(h);
87    for (i=0; i< count;i++ )
88        trade[i] = buf[i].getValue();
89
90    ierr = MPI_Send(trade, h, datatype, dest, tag, comm);
91    free(trade);
92
93    put_op(send_data);
94    ADOLC_PUT_LOCINT(count);
95    for(i=0; i < count; i++)
96       ADOLC_PUT_LOCINT(buf[i].loc());
97    ADOLC_PUT_LOCINT(count);
98    ADOLC_PUT_LOCINT(dest);
99    ADOLC_PUT_LOCINT(tag);
100
101    return ierr;
102}
103
104int ADOLC_MPI_Recv( adouble *buf,
105                   int count,
106                   ADOLC_MPI_Datatype datatype,
107                   int source, int tag,
108                   ADOLC_MPI_Comm comm
109) {
110    int i,h=count;
111    double *trade;
112    int ierr =0;
113
114
115    MPI_Status status;
116    trade = (double*) myalloc1(h);
117    ierr = MPI_Recv(trade,h, datatype, source, tag, comm, &status);
118
119    if (buf==NULL)
120       buf = new adouble[count];
121    for (i=0; i< count;i++)
122        buf[i].setValue(trade[i]);
123
124    free(trade);
125
126    put_op(receive_data);
127    ADOLC_PUT_LOCINT(count);
128    for (i=0; i< count;i++)
129       ADOLC_PUT_LOCINT(buf[i].loc());
130    ADOLC_PUT_LOCINT(count);
131    ADOLC_PUT_LOCINT(source);
132    ADOLC_PUT_LOCINT(tag);
133
134    return ierr;
135}
136int ADOLC_MPI_Bcast( adouble *buf,
137                     int count,
138                     ADOLC_MPI_Datatype datatype,
139                     int root,
140                     ADOLC_MPI_Comm comm )
141
142{
143    int i,id, ierr=0;
144    double *trade;
145
146    MPI_Comm_rank(MPI_COMM_WORLD, &id);
147
148    trade = (double*) myalloc1(count);
149
150    if ( id == root)
151       for(i= 0; i < count; i++)
152          trade[i] = buf[i].getValue();
153
154    ierr = MPI_Bcast(trade,count,datatype,root, comm);
155
156    if ( id != root){
157       if (buf==NULL)
158          buf = new adouble[count];
159       for(i=0; i< count;i++)
160          buf[i].setValue(trade[i]);
161    }
162
163    free(trade);
164
165    put_op(broadcast);
166    ADOLC_PUT_LOCINT(count);
167    for (i=0; i< count;i++)
168       ADOLC_PUT_LOCINT(buf[i].loc());
169    ADOLC_PUT_LOCINT(count);
170    ADOLC_PUT_LOCINT(root);
171    ADOLC_PUT_LOCINT(id);
172
173    return ierr;
174}
175
176int ADOLC_MPI_Reduce(
177    adouble *send_buf, adouble *rec_buf, int count, ADOLC_MPI_Datatype type,
178    ADOLC_MPI_Op op, int root, ADOLC_MPI_Comm comm)
179{
180    int i,j,id,size, ierr=0;
181
182    MPI_Comm_rank(MPI_COMM_WORLD, &id);
183    MPI_Comm_size(MPI_COMM_WORLD, &size);
184    adouble *tmp_adoubles = NULL, tmp;
185    if( id == root )
186        tmp_adoubles = new adouble[size*count];
187
188    ierr = ADOLC_MPI_Gather(send_buf,tmp_adoubles,count,type,root,comm);
189    if ( id == root){
190       if( rec_buf == NULL)
191           rec_buf = new adouble[count];
192       switch (op) {
193               case ADOLC_MPI_MAX: for(i=0; i < count; i++ ) {
194                                       tmp = tmp_adoubles[i];
195                                       for(j=1; j< size ; j++)
196                                           condassign(tmp, tmp <= tmp_adoubles[j*count+i], tmp_adoubles[j*count+i] );
197                                       rec_buf[i] = tmp;
198                                   }
199                                   break;
200               case ADOLC_MPI_MIN: for(i=0; i < count; i++ ) {
201                                      tmp = tmp_adoubles[i];
202                                      for(j=1; j< size ; j++)
203                                          condassign(tmp, tmp >= tmp_adoubles[j*count+i], tmp_adoubles[j*count+i] );
204                                      rec_buf[i] = tmp;
205                                   }
206                                   break;
207               case ADOLC_MPI_SUM: for(i=0; i < count; i++ ) {
208                                      tmp =0.;
209                                      for(j=0; j< size ; j++)
210                                         tmp += tmp_adoubles[j*count+i];
211                                       rec_buf[i] = tmp;
212                                   }
213                                   break;
214               case ADOLC_MPI_PROD:for(i=0; i < count; i++ ) {
215                                      tmp = 1.;
216                                      for(j=0; j< size ; j++)
217                                         tmp *= tmp_adoubles[j*count+i];
218                                      rec_buf[i] = tmp;
219                                    }
220                                    break;
221               default:             printf("Operation %d not yet implemented!\n",op);
222                                    break;
223       }
224       delete[] tmp_adoubles;
225    }
226
227    return ierr;
228}
229
230int ADOLC_MPI_Gather(
231    adouble *sendbuf, adouble *recvbuf, int count, ADOLC_MPI_Datatype type, int root, MPI_Comm comm)
232{
233    int i,id,size, ierr=0;
234    double *trade_s, *trade_r;
235
236    MPI_Comm_rank(MPI_COMM_WORLD, &id);
237    MPI_Comm_size(MPI_COMM_WORLD, &size);
238
239    trade_s = (double*) myalloc1(count);
240    if (id == root)
241      trade_r = (double*) myalloc1(count*size);
242    else trade_r = NULL;
243
244    for(i= 0; i < count; i++) {
245       trade_s[i] = sendbuf[i].getValue();
246    }
247    ierr = MPI_Gather(trade_s,count,type,trade_r,count,type, root, comm);
248
249    if ( id == root){
250       if( recvbuf == NULL)
251           recvbuf = new adouble[count*size];
252       for(i=0; i< count*size;i++){
253          recvbuf[i].setValue(trade_r[i]);
254          }
255    free(trade_r);
256    }
257    free( trade_s);
258
259    put_op(gather);
260    ADOLC_PUT_LOCINT(count);
261    for(i= 0; i < count; i++)
262        ADOLC_PUT_LOCINT(sendbuf[i].loc());
263    ADOLC_PUT_LOCINT(count);
264    ADOLC_PUT_LOCINT(root);
265    ADOLC_PUT_LOCINT(id);
266    ADOLC_PUT_LOCINT(count*size);
267    if( id==root) for(i=0; i < count*size;i++)
268        ADOLC_PUT_LOCINT(recvbuf[i].loc());
269    ADOLC_PUT_LOCINT(count*size);
270    ADOLC_PUT_LOCINT(root);
271    ADOLC_PUT_LOCINT(id);
272
273    return ierr;
274}
275
276int ADOLC_MPI_Scatter(
277    adouble *sendbuf, int sendcount, adouble *recvbuf,
278    int recvcount, ADOLC_MPI_Datatype type, int root, MPI_Comm comm)
279{
280    int i,id,size, ierr=0;
281    double *trade_s, *trade_r;
282
283    MPI_Comm_rank(MPI_COMM_WORLD, &id);
284    MPI_Comm_size(MPI_COMM_WORLD, &size);
285
286    trade_r = (double*) myalloc1(recvcount);
287    if (id == root)
288      trade_s = (double*) myalloc1(sendcount*size);
289    else trade_s = NULL;
290
291    if ( id == root){
292       for(i= 0; i < sendcount*size; i++)
293          trade_s[i] = sendbuf[i].getValue();
294    }
295
296    ierr = MPI_Scatter(trade_s,sendcount,type,trade_r,recvcount,type, root, comm);
297
298    if( recvbuf == NULL)
299       recvbuf = new adouble[recvcount];
300
301    for(i=0; i< recvcount;i++)
302         recvbuf[i].setValue(trade_r[i]);
303
304    free(trade_s);
305    free(trade_r);
306
307    put_op(scatter);
308    ADOLC_PUT_LOCINT(sendcount*size);
309    ADOLC_PUT_LOCINT(root);
310    ADOLC_PUT_LOCINT(id);
311    if( id == root ) {
312      for(i=0; i< sendcount*size ;i++)
313        ADOLC_PUT_LOCINT(sendbuf[i].loc());
314    }
315    ADOLC_PUT_LOCINT(sendcount*size);
316    ADOLC_PUT_LOCINT(root);
317    ADOLC_PUT_LOCINT(id);
318    ADOLC_PUT_LOCINT(recvcount);
319    for(i=0; i< recvcount;i++)
320      ADOLC_PUT_LOCINT(recvbuf[i].loc());
321    ADOLC_PUT_LOCINT(recvcount);
322
323    return ierr;
324}
325
326
327/*********************************************************************/
328/* Algorithmic Differentation Programs                               */
329
330int function(int id, int size,short tag,int m,int n,double* argument ,double* result){
331     return function_mpi(id,size,tag,m,n,argument,result);
332}
333
334int gradient(int id,int size,short tag,int n,double *x,double *y)
335{
336     return gradient_mpi(id,size,tag,n,x,y);
337}
338
339int hessian(int id,int size,short tag,int n,double *x,double **x_pp)
340{
341     return hessian_mpi(id,size,tag,n,x,x_pp);
342}
343
344int jacobian(int id,int size,short tag,int m,int n,const double *x,double **x_pp)
345{
346return jacobian_mpi(id,size,tag,m,n,x,x_pp);
347}
348
349int vec_jac(int id,int size,short tag,int m,int n,int p,double *x,double *y, double *z)
350{
351return vec_jac_mpi(id,size,tag,m,n,p,x,y,z);
352}
353
354int jac_vec(int id,int size,short tag,int m,int n,double *x,double *y,double *z)
355{
356return jac_vec_mpi(id,size,tag,m,n,x,y,z);
357}
358
359int hess_vec(int id,int size,short tag,int n,double *x,double *y,double *z)
360{
361return hess_vec_mpi(id,size,tag,n,x,y,z);
362}
363
364int lagra_hess_vec(int id,int size,short tag,int n,int p,double *x,double *y,double *t,double *z)
365{
366return lagra_hess_vec_mpi(id,size,tag,n,p,x,y,t,z);
367}
368
369void tape_doc(int id,int size,short tag, int m,int n, double *x, double *y)
370{
371return tape_doc_mpi(id,size,tag,m,n,x,y);
372}
373
374BEGIN_C_DECLS
375
376/* C - functions                                   */
377int function_mpi(int id, int size,short tag,int m,int n,double* argument ,double* result){
378     int rc =-1;
379     if( id == 0)
380          rc = zos_forward_mpi(id,size,tag,m,n,0,argument,result);
381     else
382          rc = zos_forward_mpi(id,size,tag,0,0,0,NULL,NULL);
383     return rc;
384}
385
386int gradient_mpi(int id,int size,short tag ,int n, double* x,double* result){
387
388     int rc=-1;
389     double one =1.0;
390     if( id == 0){
391          rc = zos_forward_mpi(id,size,tag,1,n,1,x,result);
392          if(rc <0){
393               printf("Failure by computing parallel gradient, process id %d!\n",id);
394               return rc;
395          }
396          rc = fos_reverse_mpi(id,size,tag,1,n,&one,result);
397     } else {
398          rc = zos_forward_mpi(id,size,tag,0,0,1,NULL,NULL);
399          if(rc <0){
400               printf("Failure by computing parallel gradient, process id %d!\n",id);
401               return rc;
402          }
403          rc = fos_reverse_mpi(id,size,tag,0,0,&one,NULL);
404     }
405     return rc;
406}
407
408int jacobian_mpi(int id, int size ,short tag ,int m,int n,const double* x,double** jacobian){
409    int rc=-1;
410    double *result = NULL , **I = NULL;
411
412    if(id == 0){
413        result = myalloc1(m);
414        if (n/2 < m) {
415            I = myallocI2(n);
416            rc = fov_forward_mpi(id,size,tag,m,n,n,x,I,result,jacobian);
417            myfreeI2(n, I);
418        } else {
419            I = myallocI2(m);
420            rc = zos_forward_mpi(id,size,tag,m,n,1,x,result);
421            if(rc <0){
422               printf("Failure by computing parallel jacobian, process id %d!\n",id);
423               return rc;
424            }
425            rc = fov_reverse_mpi(id,size,tag,m,n,m,I,jacobian);
426            myfreeI2(m, I);
427        }
428        myfree1(result);
429     } else {
430        if (n/2 < m) {
431            rc = fov_forward_mpi(id,size,tag,0,0,n,NULL,NULL,NULL,NULL);
432        } else {
433            rc = zos_forward_mpi(id,size,tag,0,0,1,NULL,NULL);
434            if(rc <0){
435               printf("Failure by computing parallel jacobian, process id %d!\n",id);
436               return rc;
437            }
438            rc = fov_reverse_mpi(id,size,tag,0,0,m,NULL,NULL);
439           }
440     }
441     return rc;
442}
443
444int hessian_mpi(int id,int size,short tag ,int n,double* x ,double** result){
445     int rc =-3,i,j;
446     if ( id == 0){
447        double *v = myalloc1(n);
448        double *w = myalloc1(n);
449        for(i=0;i<n;i++) v[i] = 0;
450        for(i=0;i<n;i++) {
451           v[i] = 1;
452           rc = hess_vec_mpi(id,size,tag, n, x, v, w);
453           if(rc <0){
454             printf("Failure by computing parallel hessian, process id %d!\n",id);
455             free(v);
456             free(w);
457             return rc;
458           }
459           for(j=0;j<=i;j++)
460            result[i][j] = w[j];
461           v[i] = 0;
462        }
463        free(v);
464        free(w);
465     } else {
466        for (i=0; i<n; i++){
467            rc = fos_forward_mpi(id,size,tag, 0,0,2,NULL,NULL,NULL,NULL);
468            if (rc <0){
469               printf("Failure by computing parallel hessian, process id %d!\n",id);
470               return rc;
471            }
472            rc = hos_reverse_mpi(id,size,tag,0,0,1, NULL,NULL);
473        }
474     }
475     return rc;
476}
477
478/* vec_jac(rank,size,tag, m, n, repeat, x[n], u[m], v[n])                             */
479int vec_jac_mpi( int id,int size,short tag,int m,int n,int repeat ,double *x,double *u,double *v){
480     int rc = -3;
481     double *y = NULL;
482     if (id == 0){
483        if(!repeat) {
484            y = myalloc1(m);
485            rc = zos_forward_mpi(id,size,tag,m,n,1, x, y);
486            if (rc <0){
487               printf("Failure by computing parallel vec_jac, process id %d!\n",id);
488               return rc;
489            }
490        }
491       MINDEC(rc, fos_reverse_mpi(id,size,tag,m,n,u,v));
492       if (!repeat) myfree1(y);
493     } else{
494        if(!repeat) {
495           rc = zos_forward_mpi(id,size,tag,0,0,1,NULL,NULL);
496        if(rc < 0) return rc;
497        }
498        rc = fos_reverse_mpi(id,size,tag,0,0,NULL,NULL);
499     }
500     return rc;
501}
502
503/* jac_vec(rank,size,tag, m, n, x[n], v[n], u[m]);                                    */
504int jac_vec_mpi(int id,int size,short tag,int m,int n,double *x,double *v, double *u){
505     int rc = -3;
506     double *y = NULL;
507     if (id == 0){
508         y = myalloc1(m);
509         rc = fos_forward_mpi(id,size,tag, m, n, 0, x, v, y, u);
510         myfree1(y);
511     } else{
512         rc = fos_forward_mpi(id,size,tag, 0, 0, 0, NULL, NULL, NULL,NULL);
513     }
514     return rc;
515}
516
517/* hess_vec(rank,size,tag, n, x[n], v[n], w[n])                                       */
518int hess_vec_mpi(int id,int size,short tag,int n,double *x,double *v,double *w){
519     double one = 1.0;
520     return lagra_hess_vec_mpi(id,size,tag,1,n,x,v,&one,w);
521}
522
523/* lagra_hess_vec(tag, m, n, x[n], v[n], u[m], w[n])                        */
524int lagra_hess_vec_mpi(int id, int size, short tag,
525                   int m,
526                   int n,
527                   double *argument,
528                   double *tangent,
529                   double *lagrange,
530                   double *result) {
531    int rc=-1;
532    int i;
533    int degree = 1;
534    int keep = degree+1;
535    double **X, *y, *y_tangent;
536
537     if (id == 0 ){
538        X = myalloc2(n,2);
539        y = myalloc1(m);
540        y_tangent = myalloc1(m);
541        rc = fos_forward_mpi(id,size,tag, m, n, keep, argument, tangent, y, y_tangent);
542        if (rc <0){
543           printf("Failure by computing parallel lagra_hess_vec, process id %d!\n",id);
544           return rc;
545        }
546        MINDEC(rc, hos_reverse_mpi(id,size,tag, m, n, degree, lagrange, X));
547        for(i = 0; i < n; ++i)
548            result[i] = X[i][1];
549        myfree1(y_tangent);
550        myfree1(y);
551        myfree2(X);
552     } else {
553        rc = fos_forward_mpi(id,size,tag, 0, 0, keep, NULL, lagrange , NULL, NULL);
554        if (rc <0){
555           printf("Failure by computing parallel lagra_hess_vec, process id %d!\n",id);
556           return rc;
557        }
558        rc = hos_reverse_mpi(id,size,tag, 0, 0, degree, lagrange, NULL );
559    }
560     return rc;
561}
562
563void tape_doc_mpi( int id,int size,short tag, int m,int n, double* x, double* y){
564     if(id==0)
565        tape_doc(id+size*tag,m,n,x,y);
566     else
567        tape_doc(id+size*tag,0,0,x,y);
568}
569
570END_C_DECLS
571
572/* That's all*/
Note: See TracBrowser for help on using the repository browser.