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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3119 by gross, Fri Aug 20 08:08:27 2010 UTC revision 3120 by gross, Mon Aug 30 10:48:00 2010 UTC
# Line 18  Line 18 
18  #include "Common.h"  #include "Common.h"
19  #include "Paso.h"  #include "Paso.h"
20    
21  INLINE void Paso_BlockOps_Cpy_1(double* R, const double* V)  void Paso_BlockOps_allMV(dim_t n_block,dim_t n,double* D,index_t* pivot,double* x,double* b);
22  {  
23     R[0]=V[0];  #define Paso_BlockOps_Cpy_1(R, V) \
24  }  {\
25       *(R+0)=*(V+0);\
26  INLINE void Paso_BlockOps_Cpy_2(double* R, const double* V)  }\
27  {  
28     R[0]=V[0];  #define Paso_BlockOps_Cpy_2(R, V) \
29     R[1]=V[1];  {\
30  }     *(R+0)=*(V+0);\
31       *(R+1)=*(V+1);\
32  INLINE void Paso_BlockOps_Cpy_3(double* R, const double* V)  }\
33  {  
34     R[0]=V[0];  #define Paso_BlockOps_Cpy_3(R, V) \
35     R[1]=V[1];  {\
36     R[2]=V[2];     *(R+0)=*(V+0);\
37  }     *(R+1)=*(V+1);\
38       *(R+2)=*(V+2);\
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) );  #define Paso_BlockOps_Cpy_N(N,R,V) \
42  }  {\
43       memcpy((void*)R, (void*)V, ( (size_t) N ) * sizeof(double) );\
44  INLINE void Paso_BlockOps_MV_1(double* R, const double* MAT, const double* V)    }\
45  {  
46     register double S1=V[0];  #define Paso_BlockOps_MV_1( R, MAT, V)   \
47     register double A11=MAT[0];  { \
48     R[0]=A11 * S1;     register double S1=*V;\
49  }     register double A11=*MAT;\
50       *R=A11 * S1;\
51    }\
52  INLINE void Paso_BlockOps_MV_2(double* R, const double* MAT, const double* V)    
53  {  
54     register double S1=V[0];  #define Paso_BlockOps_MV_2(R, MAT, V) \
55     register double S2=V[1];  { \
56     register double A11=MAT[0];     register double S1=*(V+0);\
57     register double A12=MAT[2];     register double S2=*(V+1);\
58     register double A21=MAT[1];     register double A11=*(MAT+0);\
59     register double A22=MAT[3];     register double A12=*(MAT+2);\
60     R[0]  =A11 * S1 + A12 * S2;     register double A21=*(MAT+1);\
61     R[1]=A21 * S1 + A22 * S2;     register double A22=*(MAT+3);\
62  }     *(R+0)  =A11 * S1 + A12 * S2;\
63       *(R+1)=A21 * S1 + A22 * S2;\
64    }\
65  INLINE void Paso_BlockOps_MV_3(double* R, const double* MAT, const double* V)  
66  {  
67     register double S1=V[0];  #define Paso_BlockOps_MV_3(R, MAT, V) \
68     register double S2=V[1];  { \
69     register double S3=V[2];     register double S1=*(V+0);\
70     register double A11=MAT[0];     register double S2=*(V+1);\
71     register double A21=MAT[1];     register double S3=*(V+2);\
72     register double A31=MAT[2];     register double A11=*(MAT+0);\
73     register double A12=MAT[3];     register double A21=*(MAT+1);\
74     register double A22=MAT[4];     register double A31=*(MAT+2);\
75     register double A32=MAT[5];     register double A12=*(MAT+3);\
76     register double A13=MAT[6];     register double A22=*(MAT+4);\
77     register double A23=MAT[7];     register double A32=*(MAT+5);\
78     register double A33=MAT[8];     register double A13=*(MAT+6);\
79     R[0]=A11 * S1 + A12 * S2 + A13 * S3;     register double A23=*(MAT+7);\
80     R[1]=A21 * S1 + A22 * S2 + A23 * S3;     register double A33=*(MAT+8);\
81     R[2]=A31 * S1 + A32 * S2 + A33 * S3;     *(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 */  /* this does not work in place */
87  INLINE void Paso_BlockOps_MV_N(const dim_t N, double* R, const double* MAT, const double* V)  #define Paso_BlockOps_MV_N(N, R, MAT, V) \
88  {  { \
89     dim_t p,q;     dim_t p,q; \
90     for (p=0; p<N; p++ ) {                                   for (p=0; p<N; p++ ) {    \
91        register double S= MAT[p] * V[0];          register double S= *(MAT+p) * (*(V+0)); \
92        for ( q=1; q<N; q++ ) S+=MAT[p+N*q] * V[q];                  for ( q=1; q<N; q++ ) S+=*(MAT+p+N*q) * (*(V+q));  \
93        R[p]=S;                                                  *(R+p)=S;   \
94     }                                                     }  \
95  }  } \
96    
97  INLINE void Paso_BlockOps_Solve_N(const dim_t N, double* R, const double* MAT, const index_t* pivot, const double* V)  #define Paso_BlockOps_Solve_N(N, R, MAT, pivot, V) \
98  {  { \
99     Paso_setError(TYPE_ERROR, "Paso_BlockOps_Solve_N: Right now there is support block size less than 4 only");       *(pivot+0)=0;\
100  }     Paso_setError(TYPE_ERROR, "Paso_BlockOps_Solve_N: Right now there is support block size less than 4 only"); \
101    } \
102  INLINE void Paso_BlockOps_SMV_1(double* R, const double* MAT, const double* V)    
103  {  #define Paso_BlockOps_SMV_1(R, MAT, V)   \
104     register double S1=V[0];  {\
105     register double A11=MAT[0];     register double S1=*(V+0); \
106     R[0]-=A11 * S1;     register double A11=*(MAT+0);\
107  }     *(R+0)-=A11 * S1;\
108    }\
109  INLINE void Paso_BlockOps_SMV_2(double* R, const double* MAT, const double* V)  
110  {  #define Paso_BlockOps_SMV_2(R,MAT, V) \
111     register double S1=V[0];  {\
112     register double S2=V[1];     register double S1=*(V+0); \
113     register double A11=MAT[0];     register double S2=*(V+1);\
114     register double A12=MAT[2];     register double A11=*(MAT+0);\
115     register double A21=MAT[1];     register double A12=*(MAT+2);\
116     register double A22=MAT[3];     register double A21=*(MAT+1);\
117     R[0]-=A11 * S1 + A12 * S2;     register double A22=*(MAT+3);\
118     R[1]-=A21 * S1 + A22 * S2;     *(R+0)-=A11 * S1 + A12 * S2;\
119  }     *(R+1)-=A21 * S1 + A22 * S2;\
120    }\
121  INLINE void Paso_BlockOps_SMV_3(double* R, const double* MAT, const double* V)    
122  {  #define Paso_BlockOps_SMV_3(R, MAT, V)  \
123     register double S1=V[0];  {\
124     register double S2=V[1];     register double S1=*(V+0);\
125     register double S3=V[2];     register double S2=*(V+1);\
126     register double A11=MAT[0];     register double S3=*(V+2);\
127     register double A21=MAT[1];     register double A11=*(MAT+0);\
128     register double A31=MAT[2];     register double A21=*(MAT+1);\
129     register double A12=MAT[3];     register double A31=*(MAT+2);\
130     register double A22=MAT[4];     register double A12=*(MAT+3);\
131     register double A32=MAT[5];     register double A22=*(MAT+4);\
132     register double A13=MAT[6];     register double A32=*(MAT+5);\
133     register double A23=MAT[7];     register double A13=*(MAT+6);\
134     register double A33=MAT[8];     register double A23=*(MAT+7);\
135     R[0]-=A11 * S1 + A12 * S2 + A13 * S3;     register double A33=*(MAT+8);\
136     R[1]-=A21 * S1 + A22 * S2 + A23 * S3;     *(R+0)-=A11 * S1 + A12 * S2 + A13 * S3; \
137     R[2]-=A31 * S1 + A32 * S2 + A33 * S3;     *(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 */  /* this does not work in place */
 INLINE void Paso_BlockOps_SMV_N(const dim_t N, double* R, const double* MAT, const double* V)    
 {  
    dim_t p,q;  
    for (p=0; p<N; p++ ) {                                
       register double S= MAT[p] * V[0];    
       for ( q=1; q<N; q++ ) S+=MAT[p+N*q] * V[q];            
       R[p]-=S;                                            
    }                                                  
 }  
       
 #endif /* #INC_Paso_BlockOpsOPS */  
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 */

Legend:
Removed from v.3119  
changed lines
  Added in v.3120

  ViewVC Help
Powered by ViewVC 1.1.26