/[escript]/trunk/paso/src/BlockOps.h
ViewVC logotype

Contents of /trunk/paso/src/BlockOps.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3379 - (show annotations)
Wed Nov 24 04:48:49 2010 UTC (8 years ago) by gross
File MIME type: text/plain
File size: 3770 byte(s)
some clarification on lumping
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #ifndef INC_PASO_BLOCKOPS
16 #define INC_PASO_BLOCKOPS
17
18 #include "Common.h"
19 #include "Paso.h"
20
21 void Paso_BlockOps_allMV(dim_t n_block,dim_t n,double* D,index_t* pivot,double* x);
22
23
24 /* do we need this ? */
25 #define Paso_BlockOps_Cpy_1(R, V) \
26 {\
27 *(R+0)=*(V+0);\
28 }\
29
30 #define Paso_BlockOps_Cpy_2(R, V) \
31 {\
32 *(R+0)=*(V+0);\
33 *(R+1)=*(V+1);\
34 }\
35
36 #define Paso_BlockOps_Cpy_3(R, V) \
37 {\
38 *(R+0)=*(V+0);\
39 *(R+1)=*(V+1);\
40 *(R+2)=*(V+2);\
41 }\
42
43 #define Paso_BlockOps_Cpy_N(N,R,V) \
44 {\
45 memcpy((void*)R, (void*)V, ( (size_t) N ) * sizeof(double) );\
46 }\
47
48 #define Paso_BlockOps_MV_1( R, MAT, V) \
49 { \
50 register double S1=*V;\
51 register double A11=*MAT;\
52 *R=A11 * S1;\
53 }\
54
55
56 #define Paso_BlockOps_MV_2(R, MAT, V) \
57 { \
58 register double S1=*(V+0);\
59 register double S2=*(V+1);\
60 register double A11=*(MAT+0);\
61 register double A12=*(MAT+2);\
62 register double A21=*(MAT+1);\
63 register double A22=*(MAT+3);\
64 *(R+0) = A11 * S1 + A12 * S2;\
65 *(R+1) = A21 * S1 + A22 * S2;\
66 }\
67
68
69 #define Paso_BlockOps_MV_3(R, MAT, V) \
70 { \
71 register double S1=*(V+0);\
72 register double S2=*(V+1);\
73 register double S3=*(V+2);\
74 register double A11=*(MAT+0);\
75 register double A21=*(MAT+1);\
76 register double A31=*(MAT+2);\
77 register double A12=*(MAT+3);\
78 register double A22=*(MAT+4);\
79 register double A32=*(MAT+5);\
80 register double A13=*(MAT+6);\
81 register double A23=*(MAT+7);\
82 register double A33=*(MAT+8);\
83 *(R+0)=A11 * S1 + A12 * S2 + A13 * S3; \
84 *(R+1)=A21 * S1 + A22 * S2 + A23 * S3;\
85 *(R+2)=A31 * S1 + A32 * S2 + A33 * S3;\
86 }\
87
88 /* this does not work in place */
89 #define Paso_BlockOps_MV_N(N, R, MAT, V) \
90 { \
91 dim_t p,q; \
92 for (p=0; p<N; p++ ) { \
93 register double S= *(MAT+p) * (*(V+0)); \
94 for ( q=1; q<N; q++ ) S+=*(MAT+p+N*q) * (*(V+q)); \
95 *(R+p)=S; \
96 } \
97 } \
98
99 #define Paso_BlockOps_Solve_N(N, R, MAT, pivot, V) \
100 { \
101 *(pivot+0)=0;\
102 Esys_setError(TYPE_ERROR, "Paso_BlockOps_Solve_N: Right now there is support block size less than 4 only"); \
103 } \
104
105 #define Paso_BlockOps_SMV_1(R, MAT, V) \
106 {\
107 register double S1=*(V+0); \
108 register double A11=*(MAT+0);\
109 *(R+0)-=A11 * S1;\
110 }\
111
112 #define Paso_BlockOps_SMV_2(R,MAT, V) \
113 {\
114 register double S1=*(V+0); \
115 register double S2=*(V+1);\
116 register double A11=*(MAT+0);\
117 register double A12=*(MAT+2);\
118 register double A21=*(MAT+1);\
119 register double A22=*(MAT+3);\
120 *(R+0)-=A11 * S1 + A12 * S2;\
121 *(R+1)-=A21 * S1 + A22 * S2;\
122 }\
123
124 #define Paso_BlockOps_SMV_3(R, MAT, V) \
125 {\
126 register double S1=*(V+0);\
127 register double S2=*(V+1);\
128 register double S3=*(V+2);\
129 register double A11=*(MAT+0);\
130 register double A21=*(MAT+1);\
131 register double A31=*(MAT+2);\
132 register double A12=*(MAT+3);\
133 register double A22=*(MAT+4);\
134 register double A32=*(MAT+5);\
135 register double A13=*(MAT+6);\
136 register double A23=*(MAT+7);\
137 register double A33=*(MAT+8);\
138 *(R+0)-=A11 * S1 + A12 * S2 + A13 * S3; \
139 *(R+1)-=A21 * S1 + A22 * S2 + A23 * S3;\
140 *(R+2)-=A31 * S1 + A32 * S2 + A33 * S3;\
141 }\
142
143 /* this does not work in place */
144 #define Paso_BlockOps_SMV_N(N, R, MAT, V) \
145 {\
146 dim_t p,q; \
147 for (p=0; p<N; p++ ) { \
148 register double S= *(MAT+p) * *(V+0); \
149 for ( q=1; q<N; q++ ) S+=*(MAT+p+N*q) * (*(V+q)); \
150 *(R+p)-=S; \
151 }\
152 } \
153
154 #endif /* #INC_Paso_BlockOpsOPS */

  ViewVC Help
Powered by ViewVC 1.1.26