TridiagLU  1.0
Scalable, parallel solver for tridiagonal system of equations
matops.h
Go to the documentation of this file.
1 
6 #include <stdio.h>
7 
12 #define _MatrixZero_(A,N) \
13  { \
14  int arraycounter; \
15  for (arraycounter=0; arraycounter<(N)*(N); arraycounter++) *((A)+arraycounter) = 0.0; \
16  }
17 
22 #define _MatrixMultiply_(A,B,C,N) \
23  { \
24  int matopsi,matopsj,matopsk; \
25  for (matopsi=0; matopsi<(N); matopsi++) \
26  for (matopsj=0; matopsj<(N); matopsj++) { \
27  *((C)+matopsi*(N)+matopsj) = 0; \
28  for (matopsk=0; matopsk<(N); matopsk++) *((C)+matopsi*(N)+matopsj) += ((*((A)+matopsi*(N)+matopsk)) * (*((B)+matopsk*(N)+matopsj))); \
29  } \
30  }
31 
37 #define _MatrixMultiplyNonSquare_(A,B,C,NRowA,NColA,NColB) \
38  { \
39  int matopsi,matopsj,matopsk; \
40  for (matopsi=0; matopsi<(NRowA); matopsi++) \
41  for (matopsj=0; matopsj<(NColB); matopsj++) { \
42  *((C)+matopsi*(NColB)+matopsj) = 0; \
43  for (matopsk=0; matopsk<(NColA); matopsk++) *((C)+matopsi*(NColB)+matopsj) += ((*((A)+matopsi*(NColA)+matopsk)) * (*((B)+matopsk*(NColB)+matopsj))); \
44  } \
45  }
46 
51 #define _MatrixMultiplySubtract_(C,A,B,N) \
52  { \
53  int matopsi,matopsj,matopsk; \
54  for (matopsi=0; matopsi<(N); matopsi++) \
55  for (matopsj=0; matopsj<(N); matopsj++) \
56  for (matopsk=0; matopsk<(N); matopsk++) *((C)+matopsi*(N)+matopsj) -= ((*((A)+matopsi*(N)+matopsk)) * (*((B)+matopsk*(N)+matopsj))); \
57  } \
58 
59 
64 #define _MatVecMultiply_(A,x,y,N) \
65  { \
66  int matopsi,matopsj; \
67  for (matopsi=0; matopsi<(N); matopsi++) { \
68  *((y)+matopsi) = 0; \
69  for (matopsj=0; matopsj<(N); matopsj++) *((y)+matopsi) += (*((A)+matopsi*N+matopsj) * *((x)+matopsj)); \
70  } \
71  }
72 
78 #define _MatVecMultiplySubtract_(y,A,x,N) \
79  { \
80  int matopsi,matopsj; \
81  for (matopsi=0; matopsi<N; matopsi++) \
82  for (matopsj=0; matopsj<(N); matopsj++) *((y)+matopsi) -= (*((A)+matopsi*N+matopsj) * *((x)+matopsj)); \
83  }
84 
90 #define _MatrixInvert_(A,B,N) \
91  { \
92  int matopsi,matopsj,matopsk; \
93  double matopsfactor, matopssum, matopsAc[(N)*(N)]; \
94  \
95  /* make a copy of A */ \
96  for (matopsi=0; matopsi<(N)*(N); matopsi++) *(matopsAc+matopsi) = *((A)+matopsi); \
97  \
98  /* set B as the identity matrix */ \
99  for (matopsi=0; matopsi<(N)*(N); matopsi++) *((B)+matopsi) = 0.0; \
100  for (matopsi=0; matopsi<(N) ; matopsi++) *((B)+matopsi*(N)+matopsi) = 1.0; \
101  \
102  /* LU Decomposition - Forward Sweep */ \
103  for (matopsi=0; matopsi<(N)-1; matopsi++) { \
104  if (*(matopsAc+matopsi*(N)+matopsi) == 0) fprintf(stderr,"Error in MatrixInvert(): Matrix is singular.\n"); \
105  for (matopsj=matopsi+1; matopsj<(N); matopsj++) { \
106  matopsfactor = *(matopsAc+matopsj*(N)+matopsi) / *(matopsAc+matopsi*(N)+matopsi); \
107  for (matopsk=matopsi+1; matopsk<(N); matopsk++) *(matopsAc+matopsj*(N)+matopsk) -= (matopsfactor * *(matopsAc +matopsi*(N)+matopsk)); \
108  for (matopsk=0; matopsk<matopsj; matopsk++) *((B)+matopsj*(N)+matopsk) -= (matopsfactor * *((B)+matopsi*(N)+matopsk)); \
109  } \
110  } \
111  \
112  /* LU Decomposition - Backward Sweep */ \
113  for (matopsi=(N)-1; matopsi>=0; matopsi--) { \
114  for (matopsk=0; matopsk<(N); matopsk++) { \
115  matopssum = 0.0; \
116  for (matopsj=matopsi+1; matopsj<(N); matopsj++) matopssum += (*(matopsAc+matopsi*(N)+matopsj) * *((B)+matopsj*(N)+matopsk)); \
117  *((B)+matopsi*(N)+matopsk) = (*((B)+matopsi*(N)+matopsk) - matopssum) / *(matopsAc+matopsi*(N)+matopsi); \
118  } \
119  } \
120  \
121  /* Done - B contains A^{-1} now */ \
122  }