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

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

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