Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Xxt.hpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Xxt.hpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: wrapper of functions around XXT routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #ifndef NEKTAR_LIB_UTILITIES_COMMUNICATION_XXT_HPP
37 #define NEKTAR_LIB_UTILITIES_COMMUNICATION_XXT_HPP
38 
39 #include <iostream>
40 
45 #ifdef NEKTAR_USE_MPI
47 #endif
48 
49 namespace Xxt
50 {
51 using namespace std;
52 using namespace Nektar;
53 #ifdef NEKTAR_USE_MPI
54 typedef MPI_Comm comm_ext;
55 typedef MPI_Request comm_req;
56 #else
57 typedef int comm_ext;
58 typedef int comm_req;
59 #endif
60 
61 struct comm
62 {
63  unsigned int id;
64  unsigned int np;
65  comm_ext c;
66 };
67 
69 {
70  unsigned int n, *Lrp, *Lj;
71  double *L, *D;
72 };
73 
74 struct csr_mat
75 {
76  unsigned int n, *Arp, *Aj;
77  double *A;
78 };
79 
80 struct crs_data
81 {
82 
83  /* communication */
84  struct comm comm;
85  unsigned int pcoord; /* coordinate in communication tree */
86  unsigned plevels; /* # of stages of communication */
87  signed int *pother; /* let p = pother[i], then during stage i of fan-in,
88  if p>=0, receive from p
89  if p< 0, send to (-p-1)
90  fan-out is just the reverse ...
91  on proc 0, pother is never negative
92  on others, pother is negative for the last stage only */
93  comm_req *req;
94 
95  /* separators */
96  unsigned nsep; /* number of separators */
97  unsigned int *sep_size; /* # of dofs on each separator,
98  ordered from the bottom to the top of the tree:
99  separator 0 is the bottom-most one (dofs not shared)
100  separator nsep-1 is the root of the tree */
101 
102  unsigned null_space;
103  double *share_weight;
104 
105  /* vector sizes */
106  unsigned int un; /* user's vector size */
107 
108  /* xxt_solve works with "condensed" vectors;
109  same dofs as in user's vectors, but no duplicates and no Dirichlet
110  nodes, and also ordered topologically (children before parents)
111  according to the separator tree */
112 
113  unsigned int cn; /* size of condensed vectors */
114  signed int *perm_u2c; /* permutation from user vector to condensed
115  vector, p=perm_u2c[i]; xu[i] = p=-1 ? 0 : xc[p];*/
116  unsigned int ln, sn; /* xc[0 ... ln-1] are not shared
117  (ln=sep_size[0])
118  xc[ln ... ln+sn-1] are shared
119  ln+sn = cn */
120 
121  unsigned int xn; /* # of columns of x = sum_i(sep_size[i]) -
122  sep_size[0] */
123 
124  /* data */
125  struct sparse_cholesky fac_A_ll;
126  struct csr_mat A_sl;
127  unsigned int *Xp;
128  double *X; /* column i of X starts at X[Xp[i]] */
129 
130  /* execution buffers */
131  double *vl, *vc, *vx, *combuf;
132 };
133 
134 extern "C" {
135 struct crs_data *nektar_crs_setup(unsigned int n, const unsigned long *id,
136  unsigned int nz, const unsigned int *Ai,
137  const unsigned int *Aj, const double *A,
138  unsigned int null_space,
139  const struct comm *comm);
140 void nektar_crs_solve(double *x, struct crs_data *data, double *b);
141 void nektar_crs_stats(struct crs_data *data);
142 void nektar_crs_free(struct crs_data *data);
143 }
144 
145 /**
146  * @brief Initialise the matrix-solve.
147  *
148  * On each process an array of IDs for each global degree of freedom is
149  * supplied which corresponds to a unique numbering of universal degrees of
150  * freedom. Three vectors describing the matrix are also provided. The
151  * parallel matrix solve is then set up.
152  *
153  * @param pId Array of integers providing universal IDs for each
154  * global DOF on the process.
155  * @param pAi Row indices of matrix entries
156  * @param pAj Column indices of matrix entries
157  * @param pAr Values of matrix entries
158  * @param pComm Communication object used for inter-process
159  * communication.
160  * @returns crs_data structure
161  */
162 static inline struct crs_data *Init(
163  unsigned int pRank, const Nektar::Array<OneD, unsigned long> pId,
167  const LibUtilities::CommSharedPtr &pComm)
168 {
169 #ifdef NEKTAR_USE_MPI
170  unsigned int nz = pAr.num_elements();
172  boost::dynamic_pointer_cast<LibUtilities::CommMpi>(pComm);
173  ASSERTL1(vCommMpi, "Failed to cast MPI Comm object.");
174  comm vComm;
175  MPI_Comm_dup(vCommMpi->GetComm(), &vComm.c);
176  vComm.id = vCommMpi->GetRank();
177  vComm.np = vCommMpi->GetSize();
178  crs_data *result = nektar_crs_setup(pRank, &pId[0], nz, &pAi[0], &pAj[0],
179  &pAr[0], 0, &vComm);
180  MPI_Comm_free(&vComm.c);
181  return result;
182 #else
183  return 0;
184 #endif
185 }
186 
187 /**
188  * @brief Solve the matrix system for a given input vector b.
189  */
190 static inline void Solve(Nektar::Array<OneD, NekDouble> pX,
191  struct crs_data *pCrs,
193 {
194 #ifdef NEKTAR_USE_MPI
195  if (!pCrs)
196  {
197  return;
198  }
199  nektar_crs_solve(&pX[0], pCrs, &pB[0]);
200 #endif
201 }
202 
203 /**
204  * @brief Deallocates the crs mapping data.
205  */
206 static inline void Finalise(crs_data *pCrs)
207 {
208 #ifdef NEKTAR_USE_MPI
209  int finalized;
210  MPI_Finalized(&finalized);
211  if (pCrs && !finalized)
212  {
213  nektar_crs_free(pCrs);
214  }
215 #endif
216 }
217 }
218 
219 #endif
int comm_req
Definition: Xxt.hpp:58
int comm_ext
Definition: Xxt.hpp:57
comm_req * req
Definition: Xxt.hpp:93
comm_ext c
Definition: Xxt.hpp:65
unsigned int * sep_size
Definition: Xxt.hpp:97
signed int * perm_u2c
Definition: Xxt.hpp:114
unsigned int un
Definition: Xxt.hpp:106
STL namespace.
signed int * pother
Definition: Xxt.hpp:87
static void Solve(Nektar::Array< OneD, NekDouble > pX, struct crs_data *pCrs, Nektar::Array< OneD, NekDouble > pB)
Solve the matrix system for a given input vector b.
Definition: Xxt.hpp:190
unsigned int n
Definition: Xxt.hpp:70
unsigned null_space
Definition: Xxt.hpp:102
void nektar_crs_free(struct crs_data *data)
double * A
Definition: Xxt.hpp:77
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
unsigned int xn
Definition: Xxt.hpp:121
unsigned int sn
Definition: Xxt.hpp:116
boost::shared_ptr< CommMpi > CommMpiSharedPtr
Pointer to a Communicator object.
Definition: CommMpi.h:55
double * vx
Definition: Xxt.hpp:131
Definition: Xxt.hpp:49
double * L
Definition: Xxt.hpp:71
unsigned int n
Definition: Xxt.hpp:76
unsigned nsep
Definition: Xxt.hpp:96
static struct crs_data * Init(unsigned int pRank, const Nektar::Array< OneD, unsigned long > pId, const Nektar::Array< OneD, unsigned int > pAi, const Nektar::Array< OneD, unsigned int > pAj, const Nektar::Array< OneD, NekDouble > pAr, const LibUtilities::CommSharedPtr &pComm)
Initialise the matrix-solve.
Definition: Xxt.hpp:162
A global linear system.
Definition: CommMpi.h:61
struct crs_data * nektar_crs_setup(unsigned int n, const unsigned long *id, unsigned int nz, const unsigned int *Ai, const unsigned int *Aj, const double *A, unsigned int null_space, const struct comm *comm)
unsigned int cn
Definition: Xxt.hpp:113
double * share_weight
Definition: Xxt.hpp:103
static void Finalise(crs_data *pCrs)
Deallocates the crs mapping data.
Definition: Xxt.hpp:206
unsigned plevels
Definition: Xxt.hpp:86
unsigned int * Xp
Definition: Xxt.hpp:127
unsigned int id
Definition: Xxt.hpp:63
unsigned int np
Definition: Xxt.hpp:64
double * X
Definition: Xxt.hpp:128
unsigned int pcoord
Definition: Xxt.hpp:85
void nektar_crs_stats(struct crs_data *data)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void nektar_crs_solve(double *x, struct crs_data *data, double *b)