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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26