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