paso/src/Solver_preconditioner.c


/*******************************************************
*
* Copyright (c) 2003-2009 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


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

/* Paso: SystemMatrix: sets-up the preconditioner           */

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

/* Copyrights by ACcESS Australia 2003/04 */
/* Author: Lutz Gross, l.gross@uq.edu.au */

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

#include "Paso.h"
#include "SystemMatrix.h"
#include "Solver.h"

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

/*  free space */

void Paso_Preconditioner_free(Paso_Solver_Preconditioner* in) {
    if (in!=NULL) {
      Paso_Solver_ILU_free(in->ilu);
      Paso_Solver_RILU_free(in->rilu);
      Paso_Solver_Jacobi_free(in->jacobi);
      Paso_Solver_GS_free(in->gs);
      Paso_Solver_AMG_free(in->amg);
      Paso_Solver_AMG_System_free(in->amgSystem);
      MEMFREE(in);
    }
}
/*  call the iterative solver: */

void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {

    Paso_Solver_Preconditioner* prec=NULL;
    if (A->solver==NULL) {
        /* allocate structure to hold preconditioner */
        prec=MEMALLOC(1,Paso_Solver_Preconditioner);
        prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
        if (Paso_checkPtr(prec)) return;
        prec->type=UNKNOWN;
        prec->rilu=NULL;
        prec->ilu=NULL;
        prec->jacobi=NULL;
        prec->gs=NULL;
        prec->amg=NULL;
        prec->amgSystem->amgblock1=NULL;
        prec->amgSystem->amgblock2=NULL;
        prec->amgSystem->amgblock3=NULL;
        prec->amgSystem->block1=NULL;
        prec->amgSystem->block2=NULL;
        prec->amgSystem->block3=NULL;
        A->solver=prec;
        switch (options->preconditioner) {
           default:
           case PASO_JACOBI:
              if (options->verbose) printf("Jacobi preconditioner is used.\n");
              prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
              prec->type=PASO_JACOBI;
              break;
           case PASO_ILU0:
              if (options->verbose) printf("ILU preconditioner is used.\n");
              prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
              prec->type=PASO_ILU0;
              break;
           case PASO_RILU:
              if (options->verbose) printf("RILU preconditioner is used.\n");
              prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
              prec->type=PASO_RILU;
              break;
            case PASO_GS:
              if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
              prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose);
              prec->gs->sweeps=options->sweeps;
              prec->type=PASO_GS;
              break;
            case PASO_AMG:
              if (options->verbose) printf("AMG preconditioner is used.\n");
              
              if (A->row_block_size==1) {
                prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);  
              }
              else if (A->row_block_size==2) {
                
                prec->amgSystem->block1=Paso_SparseMatrix_getBlock(A->mainBlock,1);
                prec->amgSystem->block2=Paso_SparseMatrix_getBlock(A->mainBlock,2);

                prec->amgSystem->amgblock1=Paso_Solver_getAMG(prec->amgSystem->block1,options->level_max,options);
                prec->amgSystem->amgblock2=Paso_Solver_getAMG(prec->amgSystem->block2,options->level_max,options);
                
              }
              else if (A->row_block_size==3)
              {

                prec->amgSystem->block1=Paso_SparseMatrix_getBlock(A->mainBlock,1);
                prec->amgSystem->block2=Paso_SparseMatrix_getBlock(A->mainBlock,2);
                prec->amgSystem->block3=Paso_SparseMatrix_getBlock(A->mainBlock,3);
                
                /*
                Paso_SparseMatrix_saveMM(prec->amgSystem->block1,"amgtemp1.mat");
                Paso_SparseMatrix_saveMM(prec->amgSystem->block2,"amgtemp2.mat");
                Paso_SparseMatrix_saveMM(prec->amgSystem->block3,"amgtemp3.mat");
                */
                prec->amgSystem->amgblock1=Paso_Solver_getAMG(prec->amgSystem->block1,options->level_max,options);
                prec->amgSystem->amgblock2=Paso_Solver_getAMG(prec->amgSystem->block2,options->level_max,options);
                prec->amgSystem->amgblock3=Paso_Solver_getAMG(prec->amgSystem->block3,options->level_max,options);
              }
              prec->type=PASO_AMG;
              break;
 
        }
        if (! Paso_MPIInfo_noError(A->mpi_info ) ){
            Paso_Preconditioner_free(prec);
            A->solver=NULL;
        }
    }
}

/* applies the preconditioner */
/* has to be called within a parallel reqion */
/* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
    Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
    
    switch (prec->type) {
        default:
        case PASO_JACOBI:
           Paso_Solver_solveJacobi(prec->jacobi,x,b);
           break;
        case PASO_ILU0:
           Paso_Solver_solveILU(prec->ilu,x,b);
           break;
        case PASO_RILU:
           Paso_Solver_solveRILU(prec->rilu,x,b);
           break;
         case PASO_GS:
            /* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
              We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
              x_0=0
              x_1=P^{-1}(b)
              
              b_old=b
              
              b_new=b_old+b
              x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))
                 =P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)
              b_old=b_new
              
              b_new=b_old+b
              x_3=....=P^{-1}(b_new-AP^{-1}b_old)
              b_old=b_new
              
              So for n- steps we will use loop, but we always make sure that every value calculated only once!
            */

           Paso_Solver_solveGS(prec->gs,x,b);
           if (prec->gs->sweeps>1) {
           double *bold=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
           double *bnew=MEMALLOC(prec->gs->n*prec->gs->n_block,double);
           dim_t i;
           #pragma omp parallel for private(i) schedule(static)
           for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
           
           while(prec->gs->sweeps>1) {
               #pragma omp parallel for private(i) schedule(static)
               for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
                /* Compute the residual b=b-Ax*/
               Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
               /* Go round again*/
               Paso_Solver_solveGS(prec->gs,x,bnew);
               #pragma omp parallel for private(i) schedule(static)
               for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
               prec->gs->sweeps=prec->gs->sweeps-1;
           }
           
           MEMFREE(bold);
           MEMFREE(bnew); 
           
           }
           break;
         case PASO_AMG:
            if (A->row_block_size==1) {
                Paso_Solver_solveAMG(prec->amg,x,b);
            }
            else if (A->row_block_size==2) {
                dim_t i;
                double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
                double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
                double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
                double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
                
                #pragma omp parallel for private(i) schedule(static)
                for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
                    b1[i]=b[2*i];
                    b2[i]=b[2*i+1];
                }
                
                
                Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
                Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
                
                #pragma omp parallel for private(i) schedule(static)
                for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
                    x[2*i]=x1[i];
                    x[2*i+1]=x2[i];
                }
                
                MEMFREE(x1);
                MEMFREE(x2);
                MEMFREE(b1);
                MEMFREE(b2);

            }
            else if (A->row_block_size==3) {
                dim_t i;
                
                double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
                double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
                double *x3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
                double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
                double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
                double *b3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
                
                #pragma omp parallel for private(i) schedule(static)
                for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
                    b1[i]=b[3*i];
                    b2[i]=b[3*i+1];
                    b3[i]=b[3*i+2];
                 }
                
                Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
                Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
                Paso_Solver_solveAMG(prec->amgSystem->amgblock3,x3,b3);
                
                #pragma omp parallel for private(i) schedule(static)
                for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
                    x[3*i]=x1[i];
                    x[3*i+1]=x2[i];
                    x[3*i+2]=x3[i];
                }
                
                MEMFREE(x1);
                MEMFREE(x2);
                MEMFREE(x3);
                MEMFREE(b1);
                MEMFREE(b2);
                MEMFREE(b3);
            }
        break;
    }

}
1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2009 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	/**************************************************************/
16
17	/* Paso: SystemMatrix: sets-up the preconditioner */
18
19	/**************************************************************/
20
21	/* Copyrights by ACcESS Australia 2003/04 */
22	/* Author: Lutz Gross, l.gross@uq.edu.au */
23
24	/**************************************************************/
25
26	#include "Paso.h"
27	#include "SystemMatrix.h"
28	#include "Solver.h"
29
30	/***********************************************************************************/
31
32	/* free space */
33
34	void Paso_Preconditioner_free(Paso_Solver_Preconditioner* in) {
35	if (in!=NULL) {
36	Paso_Solver_ILU_free(in->ilu);
37	Paso_Solver_RILU_free(in->rilu);
38	Paso_Solver_Jacobi_free(in->jacobi);
39	Paso_Solver_GS_free(in->gs);
40	Paso_Solver_AMG_free(in->amg);
41	Paso_Solver_AMG_System_free(in->amgSystem);
42	MEMFREE(in);
43	}
44	}
45	/* call the iterative solver: */
46
47	void Paso_Solver_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
48
49	Paso_Solver_Preconditioner* prec=NULL;
50	if (A->solver==NULL) {
51	/* allocate structure to hold preconditioner */
52	prec=MEMALLOC(1,Paso_Solver_Preconditioner);
53	prec->amgSystem=MEMALLOC(1,Paso_Solver_AMG_System);
54	if (Paso_checkPtr(prec)) return;
55	prec->type=UNKNOWN;
56	prec->rilu=NULL;
57	prec->ilu=NULL;
58	prec->jacobi=NULL;
59	prec->gs=NULL;
60	prec->amg=NULL;
61	prec->amgSystem->amgblock1=NULL;
62	prec->amgSystem->amgblock2=NULL;
63	prec->amgSystem->amgblock3=NULL;
64	prec->amgSystem->block1=NULL;
65	prec->amgSystem->block2=NULL;
66	prec->amgSystem->block3=NULL;
67	A->solver=prec;
68	switch (options->preconditioner) {
69	default:
70	case PASO_JACOBI:
71	if (options->verbose) printf("Jacobi preconditioner is used.\n");
72	prec->jacobi=Paso_Solver_getJacobi(A->mainBlock);
73	prec->type=PASO_JACOBI;
74	break;
75	case PASO_ILU0:
76	if (options->verbose) printf("ILU preconditioner is used.\n");
77	prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
78	prec->type=PASO_ILU0;
79	break;
80	case PASO_RILU:
81	if (options->verbose) printf("RILU preconditioner is used.\n");
82	prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
83	prec->type=PASO_RILU;
84	break;
85	case PASO_GS:
86	if (options->verbose) printf("Gauss-Seidel preconditioner is used.\n");
87	prec->gs=Paso_Solver_getGS(A->mainBlock,options->verbose);
88	prec->gs->sweeps=options->sweeps;
89	prec->type=PASO_GS;
90	break;
91	case PASO_AMG:
92	if (options->verbose) printf("AMG preconditioner is used.\n");
93
94	if (A->row_block_size==1) {
95	prec->amg=Paso_Solver_getAMG(A->mainBlock,options->level_max,options);
96	}
97	else if (A->row_block_size==2) {
98
99	prec->amgSystem->block1=Paso_SparseMatrix_getBlock(A->mainBlock,1);
100	prec->amgSystem->block2=Paso_SparseMatrix_getBlock(A->mainBlock,2);
101
102	prec->amgSystem->amgblock1=Paso_Solver_getAMG(prec->amgSystem->block1,options->level_max,options);
103	prec->amgSystem->amgblock2=Paso_Solver_getAMG(prec->amgSystem->block2,options->level_max,options);
104
105	}
106	else if (A->row_block_size==3)
107	{
108
109	prec->amgSystem->block1=Paso_SparseMatrix_getBlock(A->mainBlock,1);
110	prec->amgSystem->block2=Paso_SparseMatrix_getBlock(A->mainBlock,2);
111	prec->amgSystem->block3=Paso_SparseMatrix_getBlock(A->mainBlock,3);
112
113	/*
114	Paso_SparseMatrix_saveMM(prec->amgSystem->block1,"amgtemp1.mat");
115	Paso_SparseMatrix_saveMM(prec->amgSystem->block2,"amgtemp2.mat");
116	Paso_SparseMatrix_saveMM(prec->amgSystem->block3,"amgtemp3.mat");
117	*/
118	prec->amgSystem->amgblock1=Paso_Solver_getAMG(prec->amgSystem->block1,options->level_max,options);
119	prec->amgSystem->amgblock2=Paso_Solver_getAMG(prec->amgSystem->block2,options->level_max,options);
120	prec->amgSystem->amgblock3=Paso_Solver_getAMG(prec->amgSystem->block3,options->level_max,options);
121	}
122	prec->type=PASO_AMG;
123	break;
124
125	}
126	if (! Paso_MPIInfo_noError(A->mpi_info ) ){
127	Paso_Preconditioner_free(prec);
128	A->solver=NULL;
129	}
130	}
131	}
132
133	/* applies the preconditioner */
134	/* has to be called within a parallel reqion */
135	/* barrier synchronization is performed before the evaluation to make sure that the input vector is available */
136	void Paso_Solver_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
137	Paso_Solver_Preconditioner* prec=(Paso_Solver_Preconditioner*) A->solver;
138
139	switch (prec->type) {
140	default:
141	case PASO_JACOBI:
142	Paso_Solver_solveJacobi(prec->jacobi,x,b);
143	break;
144	case PASO_ILU0:
145	Paso_Solver_solveILU(prec->ilu,x,b);
146	break;
147	case PASO_RILU:
148	Paso_Solver_solveRILU(prec->rilu,x,b);
149	break;
150	case PASO_GS:
151	/* Gauss-Seidel preconditioner P=U^{-1}DL^{-1} is used here with sweeps paramenter.
152	We want to solve x_new=x_old+P^{-1}(b-Ax_old). So for fisrt 3 we will have the following:
153	x_0=0
154	x_1=P^{-1}(b)
155
156	b_old=b
157
158	b_new=b_old+b
159	x_2=x_1+P^{-1}(b-Ax_1)=P^{-1}(b)+P^{-1}(b)-P^{-1}AP^{-1}(b))
160	=P^{-1}(2b-AP^{-1}b)=P^{-1}(b_new-AP^{-1}b_old)
161	b_old=b_new
162
163	b_new=b_old+b
164	x_3=....=P^{-1}(b_new-AP^{-1}b_old)
165	b_old=b_new
166
167	So for n- steps we will use loop, but we always make sure that every value calculated only once!
168	*/
169
170	Paso_Solver_solveGS(prec->gs,x,b);
171	if (prec->gs->sweeps>1) {
172	double bold=MEMALLOC(prec->gs->nprec->gs->n_block,double);
173	double bnew=MEMALLOC(prec->gs->nprec->gs->n_block,double);
174	dim_t i;
175	#pragma omp parallel for private(i) schedule(static)
176	for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=b[i];
177
178	while(prec->gs->sweeps>1) {
179	#pragma omp parallel for private(i) schedule(static)
180	for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bnew[i]=bold[i]+b[i];
181	/* Compute the residual b=b-Ax*/
182	Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), A, x, DBLE(1), bnew);
183	/* Go round again*/
184	Paso_Solver_solveGS(prec->gs,x,bnew);
185	#pragma omp parallel for private(i) schedule(static)
186	for (i=0;i<prec->gs->n*prec->gs->n_block;++i) bold[i]=bnew[i];
187	prec->gs->sweeps=prec->gs->sweeps-1;
188	}
189
190	MEMFREE(bold);
191	MEMFREE(bnew);
192
193	}
194	break;
195	case PASO_AMG:
196	if (A->row_block_size==1) {
197	Paso_Solver_solveAMG(prec->amg,x,b);
198	}
199	else if (A->row_block_size==2) {
200	dim_t i;
201	double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
202	double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
203	double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
204	double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
205
206	#pragma omp parallel for private(i) schedule(static)
207	for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
208	b1[i]=b[2*i];
209	b2[i]=b[2*i+1];
210	}
211
212
213	Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
214	Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
215
216	#pragma omp parallel for private(i) schedule(static)
217	for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
218	x[2*i]=x1[i];
219	x[2*i+1]=x2[i];
220	}
221
222	MEMFREE(x1);
223	MEMFREE(x2);
224	MEMFREE(b1);
225	MEMFREE(b2);
226
227	}
228	else if (A->row_block_size==3) {
229	dim_t i;
230
231	double *x1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
232	double *x2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
233	double *x3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
234	double *b1=MEMALLOC(prec->amgSystem->amgblock1->n,double);
235	double *b2=MEMALLOC(prec->amgSystem->amgblock2->n,double);
236	double *b3=MEMALLOC(prec->amgSystem->amgblock3->n,double);
237
238	#pragma omp parallel for private(i) schedule(static)
239	for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
240	b1[i]=b[3*i];
241	b2[i]=b[3*i+1];
242	b3[i]=b[3*i+2];
243	}
244
245	Paso_Solver_solveAMG(prec->amgSystem->amgblock1,x1,b1);
246	Paso_Solver_solveAMG(prec->amgSystem->amgblock2,x2,b2);
247	Paso_Solver_solveAMG(prec->amgSystem->amgblock3,x3,b3);
248
249	#pragma omp parallel for private(i) schedule(static)
250	for (i=0;i<prec->amgSystem->amgblock1->n;i++) {
251	x[3*i]=x1[i];
252	x[3*i+1]=x2[i];
253	x[3*i+2]=x3[i];
254	}
255
256	MEMFREE(x1);
257	MEMFREE(x2);
258	MEMFREE(x3);
259	MEMFREE(b1);
260	MEMFREE(b2);
261	MEMFREE(b3);
262	}
263	break;
264	}
265
266	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision