paso/src/Preconditioner.c


/*******************************************************
*
* Copyright (c) 2003-2010 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 "PasoUtil.h"
#include "Preconditioner.h"

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

/*  free space */

void Paso_Preconditioner_free(Paso_Preconditioner* in) {
    if (in!=NULL) {
      Paso_Preconditioner_Smoother_free(in->jacobi);
      Paso_Preconditioner_Smoother_free(in->gs);
      Paso_Preconditioner_LocalAMG_free(in->localamg);
      Paso_Preconditioner_LocalSmoother_free(in->localamgsubstitute);
      /*********************************/
      Paso_Solver_ILU_free(in->ilu);
      Paso_Solver_RILU_free(in->rilu);
      /*********************************/
      
      MEMFREE(in);
    }
}

Paso_Preconditioner* Paso_Preconditioner_alloc(Paso_SystemMatrix* A,Paso_Options* options) {

    Paso_Preconditioner* prec=NULL;

    prec=MEMALLOC(1,Paso_Preconditioner);

    if (! Esys_checkPtr(prec)) {
        
        prec->type=UNKNOWN;
    
    
    prec->jacobi=NULL;
    prec->gs=NULL;
    prec->localamg=NULL;
    prec->localamgsubstitute=NULL;
    
    /*********************************/ 
        prec->rilu=NULL;
        prec->ilu=NULL;
    /*********************************/
    
    if (options->verbose && options->use_local_preconditioner) printf("Paso: Apply preconditioner locally only.\n");

        switch (options->preconditioner) {
           default:
           case PASO_JACOBI:
          if (options->verbose) printf("Paso: Jacobi(%d) preconditioner is used.\n",options->sweeps);
          prec->jacobi=Paso_Preconditioner_Smoother_alloc(A, TRUE, options->use_local_preconditioner, options->verbose);
              prec->type=PASO_JACOBI;
          prec->sweeps=options->sweeps;
              break;
       case PASO_GS:
          if (options->verbose) printf("Paso: Gauss-Seidel(%d) preconditioner is used.\n",options->sweeps);
          prec->gs=Paso_Preconditioner_Smoother_alloc(A, FALSE, options->use_local_preconditioner, options->verbose);
          prec->type=PASO_GS;
          prec->sweeps=options->sweeps;
          break;
       case PASO_AMLI:
       case PASO_AMG:
          if (options->verbose) printf("Paso: AMG preconditioner is used.\n");
          prec->localamg=Paso_Preconditioner_LocalAMG_alloc(A->mainBlock,1,options);
          prec->sweeps=options->sweeps;
          /* if NULL is returned (and no error) no AMG has been constructed because the system is too small or not big enough 
             we now use the Smoother as a preconditioner                                                                      */
          if ( (Esys_noError()) && (prec->localamg == NULL) ) {
         if (options->verbose) { 
            if (options->smoother == PASO_JACOBI) {
            printf("Paso: Jacobi(%d) preconditioner is used.\n",prec->sweeps);
            } else {
               printf("Paso: Gauss-Seidel(%d) preconditioner is used.\n",prec->sweeps);
            }
         }
         prec->localamgsubstitute=Paso_Preconditioner_LocalSmoother_alloc(A->mainBlock, (options->smoother == PASO_JACOBI), options->verbose);
          } else {
         if (options->verbose) printf("Paso: AMG preconditioner is used.\n");
          }
          prec->type=PASO_AMG;
          Esys_MPIInfo_noError(A->mpi_info);
          break;
          
       /***************************************************************************************/   
           case PASO_ILU0:
          if (options->verbose) printf("Paso: ILU preconditioner is used.\n");
              prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
              prec->type=PASO_ILU0;
          Esys_MPIInfo_noError(A->mpi_info);
              break;
           case PASO_RILU:
          if (options->verbose) printf("Paso: RILU preconditioner is used.\n");
              prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
          Esys_MPIInfo_noError(A->mpi_info);
              prec->type=PASO_RILU;
              break;
        }
    }
    if (! Esys_MPIInfo_noError(A->mpi_info ) ){
         Paso_Preconditioner_free(prec);
    return NULL;
    } else {
    return prec;
    }
}

/* 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_Preconditioner_solve(Paso_Preconditioner* prec, Paso_SystemMatrix* A,double* x,double* b){

    switch (prec->type) {
        default:
        case PASO_JACOBI:
       Paso_Preconditioner_Smoother_solve(A, prec->jacobi,x,b,prec->sweeps, FALSE);
           break;   
    case PASO_GS:
       Paso_Preconditioner_Smoother_solve(A, prec->gs,x,b,prec->sweeps, FALSE);
       break;
    case PASO_AMG:
       if (prec->localamg == NULL) {
          Paso_Preconditioner_LocalSmoother_solve(A->mainBlock, prec->localamgsubstitute,x,b,prec->sweeps, FALSE);
       } else {
          Paso_Preconditioner_LocalAMG_solve(A->mainBlock, prec->localamg,x,b);
       }
       break;
       
    /*=========================================================*/   
        case PASO_ILU0:
           Paso_Solver_solveILU(A->mainBlock, prec->ilu,x,b);
           break;
        case PASO_RILU:
       Paso_Solver_solveRILU(prec->rilu,x,b);
           break;

    }

}
1
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	/**************************************************************/
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 "PasoUtil.h"
29	#include "Preconditioner.h"
30
31	/***********************************************************************************/
32
33	/* free space */
34
35	void Paso_Preconditioner_free(Paso_Preconditioner* in) {
36	if (in!=NULL) {
37	Paso_Preconditioner_Smoother_free(in->jacobi);
38	Paso_Preconditioner_Smoother_free(in->gs);
39	Paso_Preconditioner_LocalAMG_free(in->localamg);
40	Paso_Preconditioner_LocalSmoother_free(in->localamgsubstitute);
41	/*********************************/
42	Paso_Solver_ILU_free(in->ilu);
43	Paso_Solver_RILU_free(in->rilu);
44	/*********************************/
45
46	MEMFREE(in);
47	}
48	}
49
50	Paso_Preconditioner* Paso_Preconditioner_alloc(Paso_SystemMatrix* A,Paso_Options* options) {
51
52	Paso_Preconditioner* prec=NULL;
53
54	prec=MEMALLOC(1,Paso_Preconditioner);
55
56	if (! Esys_checkPtr(prec)) {
57
58	prec->type=UNKNOWN;
59
60
61	prec->jacobi=NULL;
62	prec->gs=NULL;
63	prec->localamg=NULL;
64	prec->localamgsubstitute=NULL;
65
66	/*********************************/
67	prec->rilu=NULL;
68	prec->ilu=NULL;
69	/*********************************/
70
71	if (options->verbose && options->use_local_preconditioner) printf("Paso: Apply preconditioner locally only.\n");
72
73	switch (options->preconditioner) {
74	default:
75	case PASO_JACOBI:
76	if (options->verbose) printf("Paso: Jacobi(%d) preconditioner is used.\n",options->sweeps);
77	prec->jacobi=Paso_Preconditioner_Smoother_alloc(A, TRUE, options->use_local_preconditioner, options->verbose);
78	prec->type=PASO_JACOBI;
79	prec->sweeps=options->sweeps;
80	break;
81	case PASO_GS:
82	if (options->verbose) printf("Paso: Gauss-Seidel(%d) preconditioner is used.\n",options->sweeps);
83	prec->gs=Paso_Preconditioner_Smoother_alloc(A, FALSE, options->use_local_preconditioner, options->verbose);
84	prec->type=PASO_GS;
85	prec->sweeps=options->sweeps;
86	break;
87	case PASO_AMLI:
88	case PASO_AMG:
89	if (options->verbose) printf("Paso: AMG preconditioner is used.\n");
90	prec->localamg=Paso_Preconditioner_LocalAMG_alloc(A->mainBlock,1,options);
91	prec->sweeps=options->sweeps;
92	/* if NULL is returned (and no error) no AMG has been constructed because the system is too small or not big enough
93	we now use the Smoother as a preconditioner */
94	if ( (Esys_noError()) && (prec->localamg == NULL) ) {
95	if (options->verbose) {
96	if (options->smoother == PASO_JACOBI) {
97	printf("Paso: Jacobi(%d) preconditioner is used.\n",prec->sweeps);
98	} else {
99	printf("Paso: Gauss-Seidel(%d) preconditioner is used.\n",prec->sweeps);
100	}
101	}
102	prec->localamgsubstitute=Paso_Preconditioner_LocalSmoother_alloc(A->mainBlock, (options->smoother == PASO_JACOBI), options->verbose);
103	} else {
104	if (options->verbose) printf("Paso: AMG preconditioner is used.\n");
105	}
106	prec->type=PASO_AMG;
107	Esys_MPIInfo_noError(A->mpi_info);
108	break;
109
110	/***************************************************************************************/
111	case PASO_ILU0:
112	if (options->verbose) printf("Paso: ILU preconditioner is used.\n");
113	prec->ilu=Paso_Solver_getILU(A->mainBlock,options->verbose);
114	prec->type=PASO_ILU0;
115	Esys_MPIInfo_noError(A->mpi_info);
116	break;
117	case PASO_RILU:
118	if (options->verbose) printf("Paso: RILU preconditioner is used.\n");
119	prec->rilu=Paso_Solver_getRILU(A->mainBlock,options->verbose);
120	Esys_MPIInfo_noError(A->mpi_info);
121	prec->type=PASO_RILU;
122	break;
123	}
124	}
125	if (! Esys_MPIInfo_noError(A->mpi_info ) ){
126	Paso_Preconditioner_free(prec);
127	return NULL;
128	} else {
129	return prec;
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_Preconditioner_solve(Paso_Preconditioner* prec, Paso_SystemMatrix* A,double* x,double* b){
137
138	switch (prec->type) {
139	default:
140	case PASO_JACOBI:
141	Paso_Preconditioner_Smoother_solve(A, prec->jacobi,x,b,prec->sweeps, FALSE);
142	break;
143	case PASO_GS:
144	Paso_Preconditioner_Smoother_solve(A, prec->gs,x,b,prec->sweeps, FALSE);
145	break;
146	case PASO_AMG:
147	if (prec->localamg == NULL) {
148	Paso_Preconditioner_LocalSmoother_solve(A->mainBlock, prec->localamgsubstitute,x,b,prec->sweeps, FALSE);
149	} else {
150	Paso_Preconditioner_LocalAMG_solve(A->mainBlock, prec->localamg,x,b);
151	}
152	break;
153
154	/=========================================================/
155	case PASO_ILU0:
156	Paso_Solver_solveILU(A->mainBlock, prec->ilu,x,b);
157	break;
158	case PASO_RILU:
159	Paso_Solver_solveRILU(prec->rilu,x,b);
160	break;
161
162	}
163
164	}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision