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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3098 - (hide annotations)
Fri Aug 20 08:08:27 2010 UTC (10 years, 2 months ago) by gross
File MIME type: text/plain
File size: 4168 byte(s)
fix to the GaussSeidel
1 gross 3097
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 gross 3098 INLINE void Paso_BlockOps_MV_3(double* R, const double* MAT, const double* V)
66 gross 3097 {
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 gross 3098 INLINE void Paso_BlockOps_MV_N(const dim_t N, double* R, const double* MAT, const double* V)
86 gross 3097 {
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 gross 3098 INLINE void Paso_BlockOps_SMV_1(double* R, const double* MAT, const double* V)
101 gross 3097 {
102     register double S1=V[0];
103     register double A11=MAT[0];
104     R[0]-=A11 * S1;
105     }
106    
107 gross 3098 INLINE void Paso_BlockOps_SMV_2(double* R, const double* MAT, const double* V)
108 gross 3097 {
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 gross 3098 INLINE void Paso_BlockOps_SMV_N(const dim_t N, double* R, const double* MAT, const double* V)
140 gross 3097 {
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