source: issm/trunk/src/c/toolkits/plapack/patches/PlapackInvertMatrix.cpp@ 831

Last change on this file since 831 was 831, checked in by Eric.Larour, 16 years ago

Added scalapack and plapack code from cielo, not yet compiling
File size: 3.6 KB
Line 
1/*!\file PlapackInvertMatrix.cpp
2 * \brief invert petsc matrix using Plapack
3 */
4
5/* petsc: */
6#include "../../petsc/petscincludes.h"
7
8/*plapack: */
9#include "../plapackincludes.h"
10
11/* Some fortran routines: */
12#include "../../scalapack/FortranMapping.h"
13
14#undef __FUNCT__
15#define __FUNCT__ "PlapackInvertMatrix"
16void PlapackInvertMatrixLocalCleanup(PLA_Obj* pa,PLA_Template* ptempl,double** parrayA,
17 int** pidxnA,MPI_Comm* pcomm_2d);
18
19int PlapackInvertMatrix(Mat* A,Mat* inv_A,int status,int con){
20 /*inv_A does not yet exist, inv_A was just xmalloced, that's all*/
21
22 /*Error management*/
23 int i,j;
24
25
26 /*input*/
27 int mA,nA;
28 int local_mA,local_nA;
29 int lower_row,upper_row,range,this_range,this_lower_row;
30 MatType type;
31
32 /*Plapack: */
33 MPI_Datatype datatype;
34 MPI_Comm comm_2d;
35 PLA_Obj a=NULL;
36 PLA_Template templ;
37 double one=1.0;
38 int ierror;
39 int nb,nb_alg;
40 int nprows,npcols;
41 int initialized=0;
42
43 /*Petsc to Plapack: */
44 double *arrayA=NULL;
45 int* idxnA=NULL;
46 int d_nz,o_nz;
47
48 /*Feedback to client*/
49 int computation_status;
50
51 /*Verify that A is square*/
52 MatGetSize(*A,&mA,&nA);
53 MatGetLocalSize(*A,&local_mA,&local_nA);
54
55 /*Some dimensions checks: */
56 if (mA!=nA) throw ErrorException(__FUNCT__," trying to take the invert of a non-square matrix!");
57
58 /* Set default Plapack parameters */
59 //First find nprows*npcols=num_procs;
60 CyclicalFactorization(&nprows,&npcols,num_procs);
61 //nprows=num_procs;
62 //npcols=1;
63 ierror = 0;
64 nb = nA/num_procs;
65 if(nA - nb*num_procs) nb++; /* without cyclic distribution */
66
67 if (ierror){
68 PLA_Set_error_checking(ierror,PETSC_TRUE,PETSC_TRUE,PETSC_FALSE );
69 }
70 else {
71 PLA_Set_error_checking(ierror,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE );
72 }
73 nb_alg = 0;
74 if (nb_alg){
75 pla_Environ_set_nb_alg (PLA_OP_ALL_ALG,nb_alg);
76 }
77
78 /*Verify that plapack is not already initialized: */
79 if(PLA_Initialized(NULL)==TRUE)PLA_Finalize();
80 /* Create a 2D communicator */
81 PLA_Comm_1D_to_2D(MPI_COMM_WORLD,nprows,npcols,&comm_2d);
82
83 /*Initlialize plapack: */
84 PLA_Init(comm_2d);
85
86 templ = NULL;
87 PLA_Temp_create(nb, 0, &templ);
88
89 /* Use suggested nb_alg if it is not provided by user */
90 if (nb_alg == 0){
91 PLA_Environ_nb_alg(PLA_OP_PAN_PAN,templ,&nb_alg);
92 pla_Environ_set_nb_alg(PLA_OP_ALL_ALG,nb_alg);
93 }
94
95 /* Set the datatype */
96 datatype = MPI_DOUBLE;
97
98
99 /* Copy A into a*/
100 PLA_Matrix_create(datatype,mA,nA,templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&a);
101 PLA_Obj_set_to_zero(a);
102 /*Take array from A: use MatGetValues, because we are sure this routine works with
103 any matrix type.*/
104 MatGetOwnershipRange(*A,&lower_row,&upper_row);
105 upper_row--;
106 range=upper_row-lower_row+1;
107 arrayA=xmalloc(nA*sizeof(double));
108 idxnA=xmalloc(nA*sizeof(int));
109 for (i=0;i<nA;i++){
110 *(idxnA+i)=i;
111 }
112 PLA_API_begin();
113 PLA_Obj_API_open(a);
114 for (i=lower_row;i<=upper_row;i++){
115 MatGetValues(*A,1,&i,nA,idxnA,arrayA);
116 PLA_API_axpy_matrix_to_global(1,nA, &one,(void *)arrayA,1,a,i,0);
117 }
118 PLA_Obj_API_close(a);
119 PLA_API_end();
120
121 #ifdef _DEBUG_
122 PLA_Global_show("Matrix A",a," %lf","Done with A");
123 MatView(*A,PETSC_VIEWER_STDOUT_WORLD);
124 #endif
125
126 /*Call the plapack invert routine*/
127 PLA_General_invert(PLA_METHOD_INV,a);
128
129 /*Translate Plapack a into Petsc invA*/
130 MatGetType(*A,&type);
131 PlapackToPetsc(inv_A,local_mA,local_nA,mA,nA,type,a,templ,nprows,npcols,nb);
132
133 #ifdef _DEBUG_
134 PLA_Global_show("Inverse of A",a," %lf","Done...");
135 MatView(*inv_A,PETSC_VIEWER_STDOUT_WORLD);
136 #endif
137
138 /*Free ressources:*/
139 PLA_Obj_free(&a);
140 PLA_Temp_free(&templ);
141 xfree((void**)&arrayA);
142 xfree((void**)&idxnA);
143
144 /*Finalize PLAPACK*/
145 PLA_Finalize();
146 MPI_Comm_free(&comm_2d);
147
148}
Note: See TracBrowser for help on using the repository browser.