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 |
#define Paso_BlockOps_Cpy_1(R, V) \ |
24 |
{\ |
25 |
*(R+0)=*(V+0);\ |
26 |
}\ |
27 |
|
28 |
#define Paso_BlockOps_Cpy_2(R, V) \ |
29 |
{\ |
30 |
*(R+0)=*(V+0);\ |
31 |
*(R+1)=*(V+1);\ |
32 |
}\ |
33 |
|
34 |
#define Paso_BlockOps_Cpy_3(R, V) \ |
35 |
{\ |
36 |
*(R+0)=*(V+0);\ |
37 |
*(R+1)=*(V+1);\ |
38 |
*(R+2)=*(V+2);\ |
39 |
}\ |
40 |
|
41 |
#define Paso_BlockOps_Cpy_N(N,R,V) \ |
42 |
{\ |
43 |
memcpy((void*)R, (void*)V, ( (size_t) N ) * sizeof(double) );\ |
44 |
}\ |
45 |
|
46 |
#define Paso_BlockOps_MV_1( R, MAT, V) \ |
47 |
{ \ |
48 |
register double S1=*V;\ |
49 |
register double A11=*MAT;\ |
50 |
*R=A11 * S1;\ |
51 |
}\ |
52 |
|
53 |
|
54 |
#define Paso_BlockOps_MV_2(R, MAT, V) \ |
55 |
{ \ |
56 |
register double S1=*(V+0);\ |
57 |
register double S2=*(V+1);\ |
58 |
register double A11=*(MAT+0);\ |
59 |
register double A12=*(MAT+2);\ |
60 |
register double A21=*(MAT+1);\ |
61 |
register double A22=*(MAT+3);\ |
62 |
*(R+0) = A11 * S1 + A12 * S2;\ |
63 |
*(R+1) = A21 * S1 + A22 * S2;\ |
64 |
}\ |
65 |
|
66 |
|
67 |
#define Paso_BlockOps_MV_3(R, MAT, V) \ |
68 |
{ \ |
69 |
register double S1=*(V+0);\ |
70 |
register double S2=*(V+1);\ |
71 |
register double S3=*(V+2);\ |
72 |
register double A11=*(MAT+0);\ |
73 |
register double A21=*(MAT+1);\ |
74 |
register double A31=*(MAT+2);\ |
75 |
register double A12=*(MAT+3);\ |
76 |
register double A22=*(MAT+4);\ |
77 |
register double A32=*(MAT+5);\ |
78 |
register double A13=*(MAT+6);\ |
79 |
register double A23=*(MAT+7);\ |
80 |
register double A33=*(MAT+8);\ |
81 |
*(R+0)=A11 * S1 + A12 * S2 + A13 * S3; \ |
82 |
*(R+1)=A21 * S1 + A22 * S2 + A23 * S3;\ |
83 |
*(R+2)=A31 * S1 + A32 * S2 + A33 * S3;\ |
84 |
}\ |
85 |
|
86 |
/* this does not work in place */ |
87 |
#define Paso_BlockOps_MV_N(N, R, MAT, V) \ |
88 |
{ \ |
89 |
dim_t p,q; \ |
90 |
for (p=0; p<N; p++ ) { \ |
91 |
register double S= *(MAT+p) * (*(V+0)); \ |
92 |
for ( q=1; q<N; q++ ) S+=*(MAT+p+N*q) * (*(V+q)); \ |
93 |
*(R+p)=S; \ |
94 |
} \ |
95 |
} \ |
96 |
|
97 |
#define Paso_BlockOps_Solve_N(N, R, MAT, pivot, V) \ |
98 |
{ \ |
99 |
*(pivot+0)=0;\ |
100 |
Esys_setError(TYPE_ERROR, "Paso_BlockOps_Solve_N: Right now there is support block size less than 4 only"); \ |
101 |
} \ |
102 |
|
103 |
#define Paso_BlockOps_SMV_1(R, MAT, V) \ |
104 |
{\ |
105 |
register double S1=*(V+0); \ |
106 |
register double A11=*(MAT+0);\ |
107 |
*(R+0)-=A11 * S1;\ |
108 |
}\ |
109 |
|
110 |
#define Paso_BlockOps_SMV_2(R,MAT, V) \ |
111 |
{\ |
112 |
register double S1=*(V+0); \ |
113 |
register double S2=*(V+1);\ |
114 |
register double A11=*(MAT+0);\ |
115 |
register double A12=*(MAT+2);\ |
116 |
register double A21=*(MAT+1);\ |
117 |
register double A22=*(MAT+3);\ |
118 |
*(R+0)-=A11 * S1 + A12 * S2;\ |
119 |
*(R+1)-=A21 * S1 + A22 * S2;\ |
120 |
}\ |
121 |
|
122 |
#define Paso_BlockOps_SMV_3(R, MAT, V) \ |
123 |
{\ |
124 |
register double S1=*(V+0);\ |
125 |
register double S2=*(V+1);\ |
126 |
register double S3=*(V+2);\ |
127 |
register double A11=*(MAT+0);\ |
128 |
register double A21=*(MAT+1);\ |
129 |
register double A31=*(MAT+2);\ |
130 |
register double A12=*(MAT+3);\ |
131 |
register double A22=*(MAT+4);\ |
132 |
register double A32=*(MAT+5);\ |
133 |
register double A13=*(MAT+6);\ |
134 |
register double A23=*(MAT+7);\ |
135 |
register double A33=*(MAT+8);\ |
136 |
*(R+0)-=A11 * S1 + A12 * S2 + A13 * S3; \ |
137 |
*(R+1)-=A21 * S1 + A22 * S2 + A23 * S3;\ |
138 |
*(R+2)-=A31 * S1 + A32 * S2 + A33 * S3;\ |
139 |
}\ |
140 |
|
141 |
/* this does not work in place */ |
142 |
#define Paso_BlockOps_SMV_N(N, R, MAT, V) \ |
143 |
{\ |
144 |
dim_t p,q; \ |
145 |
for (p=0; p<N; p++ ) { \ |
146 |
register double S= *(MAT+p) * *(V+0); \ |
147 |
for ( q=1; q<N; q++ ) S+=*(MAT+p+N*q) * (*(V+q)); \ |
148 |
*(R+p)-=S; \ |
149 |
}\ |
150 |
} \ |
151 |
|
152 |
#endif /* #INC_Paso_BlockOpsOPS */ |