/[escript]/trunk/paso/src/Solver_applyBlockDiagonalMatrix.c
ViewVC logotype

Annotation of /trunk/paso/src/Solver_applyBlockDiagonalMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 5 months ago) by ksteube
File MIME type: text/plain
File size: 2485 byte(s)
Copyright updated in all files

1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 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 dhawcroft 631
14 ksteube 1811
15 jgs 150 /**************************************************************/
16    
17     /* Paso: apply block diagonal matrix D: x=D*b */
18    
19     /* should be called within a parallel region */
20     /* barrier synconization should be performed to make sure that the input vector available */
21    
22     /**************************************************************/
23    
24     /* Copyrights by ACcESS Australia 2003, 2004, 2005 */
25     /* Author: gross@access.edu.au */
26    
27     /**************************************************************/
28    
29 gross 700 #include "Paso.h"
30 jgs 150
31     /**************************************************************/
32    
33    
34     void Paso_Solver_applyBlockDiagonalMatrix(dim_t n_block,dim_t n,double* D,index_t* pivot,double* x,double* b) {
35     dim_t i;
36 gross 495 register dim_t i3,i9;
37     register double b0,b1,b2,D00,D10,D20,D01,D11,D21,D02,D12,D22;
38    
39 jgs 150 if (n_block==1) {
40 gross 1556 #pragma omp parallel for private(i) schedule(static)
41 gross 495 for (i=0;i<n;++i) {
42 jgs 150 x[i]=D[i]*b[i];
43     }
44     } else if (n_block==2) {
45 gross 1556 #pragma omp parallel for private(i,b0,b1,D00,D10,D01,D11,i3,i9) schedule(static)
46 gross 495 for (i=0;i<n;++i) {
47     i3=2*i;
48     i9=4*i;
49     b0=b[i3];
50     b1=b[i3+1];
51     D00=D[i9 ];
52     D10=D[i9+1];
53     D01=D[i9+2];
54     D11=D[i9+3];
55     x[i3 ]=D00*b0+D01*b1;
56     x[i3+1]=D10*b0+D11*b1;
57 jgs 150 }
58     } else if (n_block==3) {
59 gross 1556 #pragma omp parallel for private(i,b0,b1,b2,D00,D10,D20,D01,D11,D21,D02,D12,D22,i3,i9) schedule(static)
60 gross 495 for (i=0;i<n;++i) {
61     i3=3*i;
62     i9=9*i;
63     b0=b[i3];
64     b1=b[i3+1];
65     b2=b[i3+2];
66     D00=D[i9 ];
67     D10=D[i9+1];
68     D20=D[i9+2];
69     D01=D[i9+3];
70     D11=D[i9+4];
71     D21=D[i9+5];
72     D02=D[i9+6];
73     D12=D[i9+7];
74     D22=D[i9+8];
75     x[i3 ]=D00*b0+D01*b1+D02*b2;
76     x[i3+1]=D10*b0+D11*b1+D12*b2;
77     x[i3+2]=D20*b0+D21*b1+D22*b2;
78 jgs 150 }
79     }
80     return;
81     }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26