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

Last change on this file since 13975 was 13975, checked in by Mathieu Morlighem, 12 years ago

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