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

Last change on this file since 3332 was 3332, checked in by Mathieu Morlighem, 16 years ago

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