src/Solvers/Solver_applyBlockDiagonalMatrix.c

/* $Id$ */


/*
********************************************************************************
*               Copyright   2006 by ACcESS MNRF                                *
*                                                                              * 
*                 http://www.access.edu.au                                     *
*           Primary Business: Queensland, Australia                            *
*     Licensed under the Open Software License version 3.0             *
*        http://www.opensource.org/licenses/osl-3.0.php                        *
********************************************************************************
*/

/**************************************************************/

/* Paso: apply block diagonal matrix D: x=D*b               */

/* should be called within a parallel region                                              */
/* barrier synconization should be performed to make sure that the input vector available */

/**************************************************************/

/* Copyrights by ACcESS Australia 2003, 2004, 2005 */
/* Author: gross@access.edu.au */

/**************************************************************/

#include "../Paso.h"

/**************************************************************/


void Paso_Solver_applyBlockDiagonalMatrix(dim_t n_block,dim_t n,double* D,index_t* pivot,double* x,double* b) {
     dim_t i;
     register dim_t i3,i9;
     register double b0,b1,b2,D00,D10,D20,D01,D11,D21,D02,D12,D22;

     if (n_block==1) {
         #pragma omp for private(i) schedule(static)
         for (i=0;i<n;++i) {
            x[i]=D[i]*b[i];
         }
     } else if (n_block==2) {
         #pragma omp for private(i,b0,b1,D00,D10,D01,D11,i3,i9) schedule(static)
         for (i=0;i<n;++i) {
            i3=2*i;
            i9=4*i;
            b0=b[i3];
            b1=b[i3+1];
            D00=D[i9  ];
            D10=D[i9+1];
            D01=D[i9+2];
            D11=D[i9+3];
            x[i3  ]=D00*b0+D01*b1;
            x[i3+1]=D10*b0+D11*b1;
         }
     } else if (n_block==3) {
         #pragma omp for private(i,b0,b1,b2,D00,D10,D20,D01,D11,D21,D02,D12,D22,i3,i9) schedule(static)
         for (i=0;i<n;++i) {
            i3=3*i;
            i9=9*i;
            b0=b[i3];
            b1=b[i3+1];
            b2=b[i3+2];
            D00=D[i9  ];
            D10=D[i9+1];
            D20=D[i9+2];
            D01=D[i9+3];
            D11=D[i9+4];
            D21=D[i9+5];
            D02=D[i9+6];
            D12=D[i9+7];
            D22=D[i9+8];
            x[i3  ]=D00*b0+D01*b1+D02*b2;
            x[i3+1]=D10*b0+D11*b1+D12*b2;
            x[i3+2]=D20*b0+D21*b1+D22*b2;
         }
     }
     return;
}

/*
 * $Log$
 * Revision 1.2  2005/09/15 03:44:40  jgs
 * Merge of development branch dev-02 back to main trunk on 2005-09-15
 *
 * Revision 1.1.2.1  2005/09/05 06:29:50  gross
 * These files have been extracted from finley to define a stand alone libray for iterative
 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
 * has not been tested yet.
 *
 *
 */
1	/* $Id$ */
2
3
4	/*
5	********************************************************************************
6	* Copyright 2006 by ACcESS MNRF *
7	* *
8	* http://www.access.edu.au *
9	* Primary Business: Queensland, Australia *
10	* Licensed under the Open Software License version 3.0 *
11	* http://www.opensource.org/licenses/osl-3.0.php *
12	********************************************************************************
13	*/
14
15	/**************************************************************/
16
17	/* Paso: apply block diagonal matrix D: x=Db /
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	#include "../Paso.h"
30
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	register dim_t i3,i9;
37	register double b0,b1,b2,D00,D10,D20,D01,D11,D21,D02,D12,D22;
38
39	if (n_block==1) {
40	#pragma omp for private(i) schedule(static)
41	for (i=0;i<n;++i) {
42	x[i]=D[i]*b[i];
43	}
44	} else if (n_block==2) {
45	#pragma omp for private(i,b0,b1,D00,D10,D01,D11,i3,i9) schedule(static)
46	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 ]=D00b0+D01b1;
56	x[i3+1]=D10b0+D11b1;
57	}
58	} else if (n_block==3) {
59	#pragma omp for private(i,b0,b1,b2,D00,D10,D20,D01,D11,D21,D02,D12,D22,i3,i9) schedule(static)
60	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 ]=D00b0+D01b1+D02*b2;
76	x[i3+1]=D10b0+D11b1+D12*b2;
77	x[i3+2]=D20b0+D21b1+D22*b2;
78	}
79	}
80	return;
81	}
82
83	/*
84	* $Log$
85	* Revision 1.2 2005/09/15 03:44:40 jgs
86	* Merge of development branch dev-02 back to main trunk on 2005-09-15
87	*
88	* Revision 1.1.2.1 2005/09/05 06:29:50 gross
89	* These files have been extracted from finley to define a stand alone libray for iterative
90	* linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
91	* has not been tested yet.
92	*
93	*
94	*/
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision