source: issm/trunk-jpl/src/c/toolkits/plapack/patches/PlapackInvertMatrix.cpp@ 13760

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

CHG: some cleanup from cppcheck
File size: 3.2 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,range;
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 int initialized=0;
39
40 /*Petsc to Plapack: */
41 double *arrayA=NULL;
42 int* idxnA=NULL;
43
44 /*Verify that A is square*/
45 MatGetSize(*A,&mA,&nA);
46 MatGetLocalSize(*A,&local_mA,&local_nA);
47
48 /*Some dimensions checks: */
49 if (mA!=nA) _error_("trying to take the invert of a non-square matrix!");
50
51 /* Set default Plapack parameters */
52 //First find nprows*npcols=num_procs;
53 CyclicalFactorization(&nprows,&npcols,num_procs);
54 //nprows=num_procs;
55 //npcols=1;
56 ierror = 0;
57 nb = nA/num_procs;
58 if(nA - nb*num_procs) nb++; /* without cyclic distribution */
59
60 if (ierror){
61 PLA_Set_error_checking(ierror,PETSC_TRUE,PETSC_TRUE,PETSC_FALSE );
62 }
63 else {
64 PLA_Set_error_checking(ierror,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE );
65 }
66 nb_alg = 0;
67 if (nb_alg){
68 pla_Environ_set_nb_alg (PLA_OP_ALL_ALG,nb_alg);
69 }
70
71 /*Verify that plapack is not already initialized: */
72 if(PLA_Initialized(NULL)==TRUE)PLA_Finalize();
73 /* Create a 2D communicator */
74 PLA_Comm_1D_to_2D(comm,nprows,npcols,&comm_2d);
75
76 /*Initlialize plapack: */
77 PLA_Init(comm_2d);
78
79 templ = NULL;
80 PLA_Temp_create(nb, 0, &templ);
81
82 /* Use suggested nb_alg if it is not provided by user */
83 if (nb_alg == 0){
84 PLA_Environ_nb_alg(PLA_OP_PAN_PAN,templ,&nb_alg);
85 pla_Environ_set_nb_alg(PLA_OP_ALL_ALG,nb_alg);
86 }
87
88 /* Set the datatype */
89 datatype = MPI_DOUBLE;
90
91 /* Copy A into a*/
92 PLA_Matrix_create(datatype,mA,nA,templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&a);
93 PLA_Obj_set_to_zero(a);
94 /*Take array from A: use MatGetValues, because we are sure this routine works with
95 any matrix type.*/
96 MatGetOwnershipRange(*A,&lower_row,&upper_row);
97 upper_row--;
98 range=upper_row-lower_row+1;
99 arrayA = xNew<double>(nA);
100 idxnA = xNew<int>(nA);
101 for (i=0;i<nA;i++){
102 *(idxnA+i)=i;
103 }
104 PLA_API_begin();
105 PLA_Obj_API_open(a);
106 for (i=lower_row;i<=upper_row;i++){
107 MatGetValues(*A,1,&i,nA,idxnA,arrayA);
108 PLA_API_axpy_matrix_to_global(1,nA, &one,(void *)arrayA,1,a,i,0);
109 }
110 PLA_Obj_API_close(a);
111 PLA_API_end();
112
113 /*Call the plapack invert routine*/
114 PLA_General_invert(PLA_METHOD_INV,a);
115
116 /*Translate Plapack a into Petsc invA*/
117 MatGetType(*A,&type);
118 PlapackToPetsc(inv_A,local_mA,local_nA,mA,nA,type,a,templ,nprows,npcols,nb);
119
120 /*Free ressources:*/
121 PLA_Obj_free(&a);
122 PLA_Temp_free(&templ);
123 xDelete<double>(arrayA);
124 xDelete<int>(idxnA);
125
126 /*Finalize PLAPACK*/
127 PLA_Finalize();
128 MPI_Comm_free(&comm_2d);
129}
Note: See TracBrowser for help on using the repository browser.