HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
matops.h
Go to the documentation of this file.
1 
6 #ifndef _MATOPS_H_
7 #define _MATOPS_H_
8 
9 #include <stdio.h>
10 
15 #define _MatrixZero_(A,N) \
16  { \
17  int arraycounter; \
18  for (arraycounter=0; arraycounter<(N)*(N); arraycounter++) *((A)+arraycounter) = 0.0; \
19  }
20 
25 #define _MatrixMultiply_(A,B,C,N) \
26  { \
27  int matopsi,matopsj,matopsk; \
28  for (matopsi=0; matopsi<(N); matopsi++) \
29  for (matopsj=0; matopsj<(N); matopsj++) { \
30  *((C)+matopsi*(N)+matopsj) = 0; \
31  for (matopsk=0; matopsk<(N); matopsk++) *((C)+matopsi*(N)+matopsj) += ((*((A)+matopsi*(N)+matopsk)) * (*((B)+matopsk*(N)+matopsj))); \
32  } \
33  }
34 
40 #define _MatrixMultiplyNonSquare_(A,B,C,NRowA,NColA,NColB) \
41  { \
42  int matopsi,matopsj,matopsk; \
43  for (matopsi=0; matopsi<(NRowA); matopsi++) \
44  for (matopsj=0; matopsj<(NColB); matopsj++) { \
45  *((C)+matopsi*(NColB)+matopsj) = 0; \
46  for (matopsk=0; matopsk<(NColA); matopsk++) *((C)+matopsi*(NColB)+matopsj) += ((*((A)+matopsi*(NColA)+matopsk)) * (*((B)+matopsk*(NColB)+matopsj))); \
47  } \
48  }
49 
54 #define _MatrixMultiplySubtract_(C,A,B,N) \
55  { \
56  int matopsi,matopsj,matopsk; \
57  for (matopsi=0; matopsi<(N); matopsi++) \
58  for (matopsj=0; matopsj<(N); matopsj++) \
59  for (matopsk=0; matopsk<(N); matopsk++) *((C)+matopsi*(N)+matopsj) -= ((*((A)+matopsi*(N)+matopsk)) * (*((B)+matopsk*(N)+matopsj))); \
60  } \
61 
62 
67 #define _MatVecMultiply_(A,x,y,N) \
68  { \
69  int matopsi,matopsj; \
70  for (matopsi=0; matopsi<(N); matopsi++) { \
71  *((y)+matopsi) = 0; \
72  for (matopsj=0; matopsj<(N); matopsj++) *((y)+matopsi) += (*((A)+matopsi*N+matopsj) * *((x)+matopsj)); \
73  } \
74  }
75 
81 #define _MatVecMultiplySubtract_(y,A,x,N) \
82  { \
83  int matopsi,matopsj; \
84  for (matopsi=0; matopsi<N; matopsi++) \
85  for (matopsj=0; matopsj<(N); matopsj++) *((y)+matopsi) -= (*((A)+matopsi*N+matopsj) * *((x)+matopsj)); \
86  }
87 
93 #define _MatrixInvert_(A,B,N) \
94  { \
95  int matopsi,matopsj,matopsk; \
96  double matopsfactor, matopssum, matopsAc[(N)*(N)]; \
97  \
98  /* make a copy of A */ \
99  for (matopsi=0; matopsi<(N)*(N); matopsi++) *(matopsAc+matopsi) = *((A)+matopsi); \
100  \
101  /* set B as the identity matrix */ \
102  for (matopsi=0; matopsi<(N)*(N); matopsi++) *((B)+matopsi) = 0.0; \
103  for (matopsi=0; matopsi<(N) ; matopsi++) *((B)+matopsi*(N)+matopsi) = 1.0; \
104  \
105  /* LU Decomposition - Forward Sweep */ \
106  for (matopsi=0; matopsi<(N)-1; matopsi++) { \
107  if (*(matopsAc+matopsi*(N)+matopsi) == 0) fprintf(stderr,"Error in MatrixInvert(): Matrix is singular.\n"); \
108  for (matopsj=matopsi+1; matopsj<(N); matopsj++) { \
109  matopsfactor = *(matopsAc+matopsj*(N)+matopsi) / *(matopsAc+matopsi*(N)+matopsi); \
110  for (matopsk=matopsi+1; matopsk<(N); matopsk++) *(matopsAc+matopsj*(N)+matopsk) -= (matopsfactor * *(matopsAc +matopsi*(N)+matopsk)); \
111  for (matopsk=0; matopsk<matopsj; matopsk++) *((B)+matopsj*(N)+matopsk) -= (matopsfactor * *((B)+matopsi*(N)+matopsk)); \
112  } \
113  } \
114  \
115  /* LU Decomposition - Backward Sweep */ \
116  for (matopsi=(N)-1; matopsi>=0; matopsi--) { \
117  for (matopsk=0; matopsk<(N); matopsk++) { \
118  matopssum = 0.0; \
119  for (matopsj=matopsi+1; matopsj<(N); matopsj++) matopssum += (*(matopsAc+matopsi*(N)+matopsj) * *((B)+matopsj*(N)+matopsk)); \
120  *((B)+matopsi*(N)+matopsk) = (*((B)+matopsi*(N)+matopsk) - matopssum) / *(matopsAc+matopsi*(N)+matopsi); \
121  } \
122  } \
123  \
124  /* Done - B contains A^{-1} now */ \
125  }
126 
127 #endif