Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
50 namespace Xxt
51 {
52  using namespace std;
53  using namespace Nektar;
54 #ifdef NEKTAR_USE_MPI
55  typedef MPI_Comm comm_ext;
56  typedef MPI_Request comm_req;
57 #else
58  typedef int comm_ext;
59  typedef int comm_req;
60 #endif
61 
62  struct comm {
63  unsigned int id;
64  unsigned int np;
65  comm_ext c;
66  };
67 
68  struct sparse_cholesky {
69  unsigned int n, *Lrp, *Lj;
70  double *L, *D;
71  };
72 
73  struct csr_mat {
74  unsigned int n, *Arp, *Aj; double *A;
75  };
76 
77  struct crs_data {
78 
79  /* communication */
80  struct comm comm;
81  unsigned int pcoord;/* coordinate in communication tree */
82  unsigned plevels; /* # of stages of communication */
83  signed int *pother; /* let p = pother[i], then during stage i of fan-in,
84  if p>=0, receive from p
85  if p< 0, send to (-p-1)
86  fan-out is just the reverse ...
87  on proc 0, pother is never negative
88  on others, pother is negative for the last stage only */
89  comm_req *req;
90 
91  /* separators */
92  unsigned nsep; /* number of separators */
93  unsigned int *sep_size; /* # of dofs on each separator,
94  ordered from the bottom to the top of the tree:
95  separator 0 is the bottom-most one (dofs not shared)
96  separator nsep-1 is the root of the tree */
97 
98  unsigned null_space;
99  double *share_weight;
100 
101  /* vector sizes */
102  unsigned int un; /* user's vector size */
103 
104  /* xxt_solve works with "condensed" vectors;
105  same dofs as in user's vectors, but no duplicates and no Dirichlet
106  nodes, and also ordered topologically (children before parents)
107  according to the separator tree */
108 
109  unsigned int cn; /* size of condensed vectors */
110  signed int *perm_u2c; /* permutation from user vector to condensed
111  vector, p=perm_u2c[i]; xu[i] = p=-1 ? 0 : xc[p];*/
112  unsigned int ln, sn; /* xc[0 ... ln-1] are not shared
113  (ln=sep_size[0])
114  xc[ln ... ln+sn-1] are shared
115  ln+sn = cn */
116 
117  unsigned int xn; /* # of columns of x = sum_i(sep_size[i]) -
118  sep_size[0] */
119 
120  /* data */
121  struct sparse_cholesky fac_A_ll;
122  struct csr_mat A_sl;
123  unsigned int *Xp; double *X; /* column i of X starts at X[Xp[i]] */
124 
125  /* execution buffers */
126  double *vl, *vc, *vx, *combuf;
127  };
128 
129 
130  extern "C"
131  {
132  struct crs_data *nektar_crs_setup(
133  unsigned int n, const unsigned long *id,
134  unsigned int nz, const unsigned int *Ai, const unsigned int *Aj,
135  const double *A, unsigned int null_space, const struct comm *comm);
136  void nektar_crs_solve(double *x, struct crs_data *data, double *b);
137  void nektar_crs_stats(struct crs_data *data);
138  void nektar_crs_free(struct crs_data *data);
139  }
140 
141  /**
142  * @brief Initialise the matrix-solve.
143  *
144  * On each process an array of IDs for each global degree of freedom is
145  * supplied which corresponds to a unique numbering of universal degrees of
146  * freedom. Three vectors describing the matrix are also provided. The
147  * parallel matrix solve is then set up.
148  *
149  * @param pId Array of integers providing universal IDs for each
150  * global DOF on the process.
151  * @param pAi Row indices of matrix entries
152  * @param pAj Column indices of matrix entries
153  * @param pAr Values of matrix entries
154  * @param pComm Communication object used for inter-process
155  * communication.
156  * @returns crs_data structure
157  */
158  static inline struct crs_data* Init ( unsigned int pRank,
163  const LibUtilities::CommSharedPtr& pComm)
164  {
165 #ifdef NEKTAR_USE_MPI
166  unsigned int nz = pAr.num_elements();
167  LibUtilities::CommMpiSharedPtr vCommMpi = boost::dynamic_pointer_cast<LibUtilities::CommMpi> (pComm);
168  ASSERTL1(vCommMpi, "Failed to cast MPI Comm object.");
169  comm vComm;
170  MPI_Comm_dup(vCommMpi->GetComm(), &vComm.c);
171  vComm.id = vCommMpi->GetRank();
172  vComm.np = vCommMpi->GetSize();
173  crs_data* result = nektar_crs_setup(pRank, &pId[0], nz, &pAi[0],
174  &pAj[0], &pAr[0], 0, &vComm);
175  MPI_Comm_free(&vComm.c);
176  return result;
177 #else
178  return 0;
179 #endif
180  }
181 
182 
183  /**
184  * @brief Solve the matrix system for a given input vector b.
185  */
186  static inline void Solve ( Nektar::Array<OneD, NekDouble> pX,
187  struct crs_data* pCrs,
189  {
190 #ifdef NEKTAR_USE_MPI
191  if (!pCrs) {
192  return;
193  }
194  nektar_crs_solve(&pX[0], pCrs, &pB[0]);
195 #endif
196  }
197 
198 
199  /**
200  * @brief Deallocates the crs mapping data.
201  */
202  static inline void Finalise (crs_data* pCrs)
203  {
204 #ifdef NEKTAR_USE_MPI
205  int finalized;
206  MPI_Finalized(&finalized);
207  if (pCrs && !finalized)
208  {
209  nektar_crs_free(pCrs);
210  }
211 #endif
212  }
213 }
214 
215 #endif
int comm_req
Definition: Xxt.hpp:59
int comm_ext
Definition: Xxt.hpp:58
comm_req * req
Definition: Xxt.hpp:89
comm_ext c
Definition: Xxt.hpp:65
unsigned int * sep_size
Definition: Xxt.hpp:93
signed int * perm_u2c
Definition: Xxt.hpp:110
unsigned int un
Definition: Xxt.hpp:102
STL namespace.
signed int * pother
Definition: Xxt.hpp:83
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:186
unsigned int n
Definition: Xxt.hpp:69
unsigned null_space
Definition: Xxt.hpp:98
void nektar_crs_free(struct crs_data *data)
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:53
unsigned int xn
Definition: Xxt.hpp:117
unsigned int sn
Definition: Xxt.hpp:112
boost::shared_ptr< CommMpi > CommMpiSharedPtr
Pointer to a Communicator object.
Definition: CommMpi.h:55
double * vx
Definition: Xxt.hpp:126
Definition: Xxt.hpp:50
double * L
Definition: Xxt.hpp:70
unsigned int n
Definition: Xxt.hpp:74
unsigned nsep
Definition: Xxt.hpp:92
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:158
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:109
double * share_weight
Definition: Xxt.hpp:99
static void Finalise(crs_data *pCrs)
Deallocates the crs mapping data.
Definition: Xxt.hpp:202
unsigned plevels
Definition: Xxt.hpp:82
unsigned int * Xp
Definition: Xxt.hpp:123
unsigned int id
Definition: Xxt.hpp:63
unsigned int np
Definition: Xxt.hpp:64
unsigned int pcoord
Definition: Xxt.hpp:81
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:191
void nektar_crs_solve(double *x, struct crs_data *data, double *b)