paso/src/Solver_AMG.c


/*******************************************************
*
* Copyright (c) 2003-2008 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: AMG preconditioner with reordering                 */

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

/* Author: artak@access.edu.au                                */

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

#include "Paso.h"
#include "Solver.h"
#include "PasoUtil.h"
#include "UMFPACK.h"
#include "Pattern_coupling.h"

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

/* free all memory used by AMG                                */

void Paso_Solver_AMG_free(Paso_Solver_AMG * in) {
     if (in!=NULL) {
        Paso_Solver_AMG_free(in->AMG_of_Schur);
        Paso_Solver_GS_free(in->GS);
        MEMFREE(in->inv_A_FF);
        MEMFREE(in->A_FF_pivot);
        Paso_SparseMatrix_free(in->A_FC);
        Paso_SparseMatrix_free(in->A_CF);
        MEMFREE(in->rows_in_F);
        MEMFREE(in->rows_in_C);
        MEMFREE(in->mask_F);
        MEMFREE(in->mask_C);
        MEMFREE(in->x_F);
        MEMFREE(in->b_F);
        MEMFREE(in->x_C);
        MEMFREE(in->b_C);
        MEMFREE(in);
     }
}

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

/*   constructs the block-block factorization of 

        [ A_FF A_FC ]
   A_p= 
        [ A_CF A_FF ]

to 

  [      I         0  ]  [ A_FF 0 ] [ I    invA_FF*A_FF ]  
  [ A_CF*invA_FF   I  ]  [   0  S ] [ 0          I      ] 


   where S=A_FF-ACF*invA_FF*A_FC within the shape of S

   then AMG is applied to S again until S becomes empty 

*/
Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,bool_t verbose,dim_t level) {
  Paso_Solver_AMG* out=NULL;
  dim_t n=A_p->numRows;
  dim_t n_block=A_p->row_block_size;
  index_t* mis_marker=NULL;  
  index_t* counter=NULL;
  index_t iPtr,*index, *where_p, iPtr_s;
  dim_t i,k,j,j0;
  Paso_SparseMatrix * schur=NULL;
  Paso_SparseMatrix * schur_withFillIn=NULL;
  double time0,time1,time2,S;
  
  /* identify independend set of rows/columns */
  mis_marker=TMPMEMALLOC(n,index_t);
  counter=TMPMEMALLOC(n,index_t);
  out=MEMALLOC(1,Paso_Solver_AMG);
  out->AMG_of_Schur=NULL;
  out->inv_A_FF=NULL;
  out->A_FF_pivot=NULL;
  out->A_FC=NULL;
  out->A_CF=NULL;
  out->rows_in_F=NULL;
  out->rows_in_C=NULL;
  out->mask_F=NULL;
  out->mask_C=NULL;
  out->x_F=NULL;
  out->b_F=NULL;
  out->x_C=NULL;
  out->b_C=NULL;
  out->A=Paso_SparseMatrix_getReference(A_p);
  out->GS=Paso_Solver_getGS(A_p,verbose);
  out->GS->sweeps=4;
  out->level=level;
  
  if ( !(Paso_checkPtr(mis_marker) || Paso_checkPtr(out) || Paso_checkPtr(counter) ) ) {
     /* identify independend set of rows/columns */
     time0=Paso_timer();
     #pragma omp parallel for private(i) schedule(static)
     for (i=0;i<n;++i) mis_marker[i]=-1;
     /*Paso_Pattern_RS(A_p,mis_marker,0.25);*/
     /*Paso_Pattern_Aggregiation(A_p,mis_marker,0.5);*/
     Paso_Pattern_coup(A_p,mis_marker,0.05);
     time2=Paso_timer()-time0;
     if (Paso_noError()) {
        #pragma omp parallel for private(i) schedule(static)
        for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
        out->n=n;
        out->n_block=n_block;
        out->n_F=Paso_Util_cumsum(n,counter);
        out->mask_F=MEMALLOC(n,index_t);
        out->rows_in_F=MEMALLOC(out->n_F,index_t);
        out->inv_A_FF=MEMALLOC(n_block*n_block*out->n_F,double);
        out->A_FF_pivot=NULL; /* later use for block size>3 */
        if (! (Paso_checkPtr(out->mask_F) || Paso_checkPtr(out->inv_A_FF) || Paso_checkPtr(out->rows_in_F) ) ) {
           #pragma omp parallel 
           {
              /* creates an index for F from mask */
              #pragma omp for private(i) schedule(static)
              for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;
              #pragma omp for private(i) schedule(static)
              for (i = 0; i < n; ++i) {
                 if  (mis_marker[i]) {
                        out->rows_in_F[counter[i]]=i;
                        out->mask_F[i]=counter[i];
                 } else {
                        out->mask_F[i]=-1;
                 }
              }
              /* Compute row-sum for getting rs(A_FF)*/
              #pragma omp for private(i,iPtr) schedule(static)
              for (i = 0; i < out->n_F; ++i) {
                out->inv_A_FF[i]=0;
                for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) {
                 out->inv_A_FF[i]+=A_p->val[iPtr];
                }
              }
              
/*          for (i=0;i<out->n_F;++i) {
          fprintf(stderr,"GET INV_A_FF %g %d of %d (%d)\n",out->inv_A_FF[i],i,out->rows_in_F[i],(out->inv_A_FF[i]==0));
          }
          fprintf(stderr,"END\n"); 
*/
              #pragma omp for private(i, where_p,iPtr,index) schedule(static)
              for (i = 0; i < out->n_F; i++) {
                /* find main diagonal */
                iPtr=A_p->pattern->ptr[out->rows_in_F[i]];
                index=&(A_p->pattern->index[iPtr]);
                where_p=(index_t*)bsearch(&out->rows_in_F[i],
                                        index,
                                        A_p->pattern->ptr[out->rows_in_F[i] + 1]-A_p->pattern->ptr[out->rows_in_F[i]],
                                        sizeof(index_t),
                                        Paso_comparIndex);
                if (where_p==NULL) {
                    Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
                } else {
                    iPtr+=(index_t)(where_p-index);
                    /* get inverse of A_FF block: */
                      S=out->inv_A_FF[i];
                      if (ABS(A_p->val[iPtr])>0.) {
                            if(ABS(S)>0.)
                               out->inv_A_FF[i]=1./S;
                            else
                            {
                               fprintf(stderr,"ROWSUM OF ROW %d->%d is 0\n",i,out->rows_in_F[i]);
                               S=0;
                               for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) {
                               if(A_p->val[iPtr]!=0) {
                                fprintf(stderr,"A_p[%d,%d]=%f ",out->rows_in_F[i],A_p->pattern->index[iPtr],A_p->val[iPtr]);
                                S+=A_p->val[iPtr];
                               }
                              }
                              fprintf(stderr,"\n SUMMMMMMM %f\n",S);
                            }
                      } else {
                            Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getAMG: Break-down in AMG decomposition: non-regular main diagonal block.");
                      }
                } 
              }
           } /* end parallel region */

           if( Paso_noError()) {
              /* if there are no nodes in the coarse level there is no more work to do */
              out->n_C=n-out->n_F;
              if (level<2) {
              /*if (out->n_C>10) {*/
                   out->rows_in_C=MEMALLOC(out->n_C,index_t);
                   out->mask_C=MEMALLOC(n,index_t);
                   if (! (Paso_checkPtr(out->mask_C) || Paso_checkPtr(out->rows_in_C) ) ) {
                       /* creates an index for C from mask */
                       #pragma omp parallel for private(i) schedule(static)
                       for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];
                       Paso_Util_cumsum(n,counter);
                       #pragma omp parallel
                       {
                          #pragma omp for private(i) schedule(static)
                          for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1;
                          #pragma omp for private(i) schedule(static)
                          for (i = 0; i < n; ++i) {
                             if  (! mis_marker[i]) {
                                out->rows_in_C[counter[i]]=i;
                                out->mask_C[i]=counter[i];
                             } else {
                                out->mask_C[i]=-1;
                             }
                          }
                      } /* end parallel region */
                      /* get A_CF block: */
                      out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);
                      if (Paso_noError()) {
                         /* get A_FC block: */
                         out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);
                         /* get A_CC block: */
                         if (Paso_noError()) {
                            schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);
                            
                            /*find the pattern of the schur complement with fill in*/
                            schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1);

                            /* copy values over*/ 
                            #pragma omp for private(i,iPtr,iPtr_s,j,j0) schedule(static)
                            for (i = 0; i < schur_withFillIn->numRows; ++i) {
                              for (iPtr=schur_withFillIn->pattern->ptr[i];iPtr<schur_withFillIn->pattern->ptr[i + 1]; ++iPtr) {
                                j=schur_withFillIn->pattern->index[iPtr];
                                schur_withFillIn->val[iPtr]=0.;
                                for (iPtr_s=schur->pattern->ptr[i];iPtr_s<schur->pattern->ptr[i + 1]; ++iPtr_s){
                                    j0=schur->pattern->index[iPtr_s];
                                    if (j==j0) {
                                      schur_withFillIn->val[iPtr]=schur->val[iPtr_s];
                                      break;
                                    }
                                }
                              }
                            }
                           time0=Paso_timer()-time0;
                           
                            if (Paso_noError()) {
                                time1=Paso_timer();
                                /* update A_CC block to get Schur complement and then apply AMG to it */
                                Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);
                                time1=Paso_timer()-time1;
                                out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level+1);
                                
                                /*Paso_SparseMatrix_free(schur);*/
                                /* Paso_SparseMatrix_free(schur_withFillIn);*/
                            }
                            /* allocate work arrays for AMG application */
                            if (Paso_noError()) {
                              out->x_F=MEMALLOC(n_block*out->n_F,double);
                              out->b_F=MEMALLOC(n_block*out->n_F,double);
                              out->x_C=MEMALLOC(n_block*out->n_C,double);
                              out->b_C=MEMALLOC(n_block*out->n_C,double);

                              if (! (Paso_checkPtr(out->x_F) || Paso_checkPtr(out->b_F) || Paso_checkPtr(out->x_C) || Paso_checkPtr(out->b_C) ) ) {
                                  #pragma omp parallel 
                                  {
                                    #pragma omp for private(i,k) schedule(static)
                                    for (i = 0; i < out->n_F; ++i) {
                                          for (k=0; k<n_block;++k) {
                                             out->x_F[i*n_block+k]=0.;
                                             out->b_F[i*n_block+k]=0.;
                                          }
                                    }
                                    #pragma omp for private(i,k) schedule(static)
                                    for (i = 0; i < out->n_C; ++i) {
                                        for (k=0; k<n_block;++k) {
                                          out->x_C[i*n_block+k]=0.;
                                          out->b_C[i*n_block+k]=0.;
                                        }
                                    }
                                  } /* end parallel region */
                              }
                            }
                         }
                     }
                 }
              }
           }
        }
     }
  }
  TMPMEMFREE(mis_marker);
  TMPMEMFREE(counter);
  if (Paso_noError()) {
      if (verbose) {
         printf("AMG: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);
         if (out->n_C>10) {
            printf("timing: AMG: MIS/reordering/elemination : %e/%e/%e\n",time2,time0,time1);
         } else {
            printf("timing: AMG: MIS: %e\n",time2); 
         }
     }
     return out;
  } else  {
     Paso_Solver_AMG_free(out);
     return NULL;
  }
}

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

/* apply AMG precondition b-> x                               

     in fact it solves 

  [      I         0  ]  [ A_FF 0 ] [ I    invA_FF*A_FF ]  [ x_F ]  = [b_F]
  [ A_CF*invA_FF   I  ]  [   0  S ] [ 0          I      ]  [ x_C ]  = [b_C]

 in the form 

   b->[b_F,b_C] 
   x_F=invA_FF*b_F
   b_C=b_C-A_CF*x_F
   x_C=AMG(b_C)
   b_F=b_F-A_FC*x_C
   x_F=invA_FF*b_F
   x<-[x_F,x_C]

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

*/

void Paso_Solver_solveAMG(Paso_Solver_AMG * amg, double * x, double * b) {
     dim_t i,k,oldsweeps,j;
     double *r=MEMALLOC(amg->n,double);
     /*Paso_Solver_GS* GS=NULL;*/
     double *bold=MEMALLOC(amg->n,double);
     double *bnew=MEMALLOC(amg->n,double);
     double *x0=MEMALLOC(amg->n,double);
     index_t iPtr;

     if (amg->level==2) {
     /*if (amg->n_C<=10) {*/
        Paso_UMFPACK1(amg->A,x,b,0);
     } else {
     
        /* presmoothing */
         Paso_Solver_solveGS(amg->GS,x,b);
         oldsweeps=amg->GS->sweeps;
         if (amg->GS->sweeps>1) {
           
           #pragma omp parallel for private(i) schedule(static)
           for (i=0;i<amg->GS->n;++i) bold[i]=b[i];
           
           while(amg->GS->sweeps>1) {
               #pragma omp parallel for private(i) schedule(static)
               for (i=0;i<amg->GS->n;++i) bnew[i]=bold[i]+b[i];
                /* Compute the residual b=b-Ax*/
               Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), amg->A, x, DBLE(1), bnew);
               /* Go round again*/
               Paso_Solver_solveGS(amg->GS,x,bnew);
               #pragma omp parallel for private(i) schedule(static)
               for (i=0;i<amg->GS->n;++i) bold[i]=bnew[i];
               amg->GS->sweeps=amg->GS->sweeps-1;
           }
           }
           amg->GS->sweeps=oldsweeps;
        /* end of presmoothing */
        
         #pragma omp parallel for private(i) schedule(static)
         for (i=0;i<amg->n;++i) r[i]=b[i];
         
         /*r=b-Ax*/ 
         Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
       
        /* b->[b_F,b_C]     */
        #pragma omp parallel for private(i) schedule(static)
        for (i=0;i<amg->n_F;++i) amg->b_F[i]=r[amg->rows_in_F[i]];
        
        #pragma omp parallel for private(i) schedule(static)
        for (i=0;i<amg->n_C;++i) amg->b_C[i]=r[amg->rows_in_C[i]];

        /* x_F=invA_FF*b_F  */
        Paso_Solver_applyBlockDiagonalMatrix(1,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F);
        
        /* b_C=b_C-A_CF*x_F */
        Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_CF,amg->x_F,1.,amg->b_C);
        
        /* x_C=AMG(b_C)     */
        Paso_Solver_solveAMG(amg->AMG_of_Schur,amg->x_C,amg->b_C);
        
        /* b_F=b_F-A_FC*x_C */
        Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_FC,amg->x_C,1.,amg->b_F);
        /* x_F=invA_FF*b_F  */
        Paso_Solver_applyBlockDiagonalMatrix(1,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F);
        /* x<-[x_F,x_C]     */
        
        #pragma omp parallel for private(i) schedule(static)
        for (i=0;i<amg->n;++i) {
            if (amg->mask_C[i]>-1) {
                 x[i]+=amg->x_C[amg->mask_C[i]];
            } else {
                 x[i]+=amg->x_F[amg->mask_F[i]];
            }
        }
        

     /*postsmoothing*/
     #pragma omp parallel for private(i) schedule(static)
     for (i=0;i<amg->n;++i) r[i]=b[i];
     
     /*for (i=0;i<amg->n;++i) fprintf(stderr,"%lf %lf %d %d %d BBB\n",r[i],x[i],amg->n,amg->A->numCols, amg->level);*/
     
     /*
     for (i = 0; i < amg->A->numRows; ++i) {
          for (iPtr=amg->A->pattern->ptr[i];iPtr<amg->A->pattern->ptr[i + 1]; ++iPtr) {
              j=amg->A->pattern->index[iPtr];
              fprintf(stderr,"A[%d,%d]=%lf ",i,j,amg->A->val[iPtr]);
          } 
          fprintf(stderr,"\n");
     }*/
     
     /*r=b-Ax*/ 
     Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
     Paso_Solver_solveGS(amg->GS,x0,r);
     
     #pragma omp parallel for private(i) schedule(static)
     for (i=0;i<amg->n;++i) {
      /*fprintf(stderr,"%lf %lf %lf\n",x[i],x0[i],x[i]+x0[i]);*/
      x[i]+=x0[i];
     }
     
     /*end of postsmoothing*/
     
     /*Paso_Solver_GS_free(amg->GS);*/
     }
     MEMFREE(r);
     MEMFREE(bold);
     MEMFREE(bnew);
     MEMFREE(x0);
    return;
}

/*
 * $Log$
 *
 */
1	jfenwick	1851
2			/*******************************************************
3			*
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
14
15			/**************************************************************/
16
17			/* Paso: AMG preconditioner with reordering */
18
19			/**************************************************************/
20
21			/* Author: artak@access.edu.au */
22
23			/**************************************************************/
24
25			#include "Paso.h"
26			#include "Solver.h"
27			#include "PasoUtil.h"
28	artak	1931	#include "UMFPACK.h"
29	phornby	1917	#include "Pattern_coupling.h"
30	jfenwick	1851
31			/**************************************************************/
32
33			/* free all memory used by AMG */
34
35			void Paso_Solver_AMG_free(Paso_Solver_AMG * in) {
36			if (in!=NULL) {
37			Paso_Solver_AMG_free(in->AMG_of_Schur);
38	artak	1939	Paso_Solver_GS_free(in->GS);
39	jfenwick	1851	MEMFREE(in->inv_A_FF);
40			MEMFREE(in->A_FF_pivot);
41			Paso_SparseMatrix_free(in->A_FC);
42			Paso_SparseMatrix_free(in->A_CF);
43			MEMFREE(in->rows_in_F);
44			MEMFREE(in->rows_in_C);
45			MEMFREE(in->mask_F);
46			MEMFREE(in->mask_C);
47			MEMFREE(in->x_F);
48			MEMFREE(in->b_F);
49			MEMFREE(in->x_C);
50			MEMFREE(in->b_C);
51			MEMFREE(in);
52			}
53			}
54
55			/**************************************************************/
56
57			/* constructs the block-block factorization of
58
59			[ A_FF A_FC ]
60			A_p=
61			[ A_CF A_FF ]
62
63			to
64
65			[ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ]
66			[ A_CF*invA_FF I ] [ 0 S ] [ 0 I ]
67
68
69			where S=A_FF-ACFinvA_FFA_FC within the shape of S
70
71			then AMG is applied to S again until S becomes empty
72
73			*/
74	artak	1890	Paso_Solver_AMG* Paso_Solver_getAMG(Paso_SparseMatrix *A_p,bool_t verbose,dim_t level) {
75	jfenwick	1851	Paso_Solver_AMG* out=NULL;
76			dim_t n=A_p->numRows;
77			dim_t n_block=A_p->row_block_size;
78			index_t* mis_marker=NULL;
79	artak	1881	index_t* counter=NULL;
80			index_t iPtr,index, where_p, iPtr_s;
81			dim_t i,k,j,j0;
82	jfenwick	1851	Paso_SparseMatrix * schur=NULL;
83	artak	1881	Paso_SparseMatrix * schur_withFillIn=NULL;
84	phornby	1917	double time0,time1,time2,S;
85	artak	1890
86	jfenwick	1851	/* identify independend set of rows/columns */
87			mis_marker=TMPMEMALLOC(n,index_t);
88			counter=TMPMEMALLOC(n,index_t);
89			out=MEMALLOC(1,Paso_Solver_AMG);
90			out->AMG_of_Schur=NULL;
91			out->inv_A_FF=NULL;
92			out->A_FF_pivot=NULL;
93			out->A_FC=NULL;
94			out->A_CF=NULL;
95			out->rows_in_F=NULL;
96			out->rows_in_C=NULL;
97			out->mask_F=NULL;
98			out->mask_C=NULL;
99			out->x_F=NULL;
100			out->b_F=NULL;
101			out->x_C=NULL;
102			out->b_C=NULL;
103	artak	1890	out->A=Paso_SparseMatrix_getReference(A_p);
104	artak	1939	out->GS=Paso_Solver_getGS(A_p,verbose);
105	artak	1954	out->GS->sweeps=4;
106	artak	1890	out->level=level;
107	artak	1881
108	jfenwick	1851	if ( !(Paso_checkPtr(mis_marker) \|\| Paso_checkPtr(out) \|\| Paso_checkPtr(counter) ) ) {
109			/* identify independend set of rows/columns */
110			time0=Paso_timer();
111			#pragma omp parallel for private(i) schedule(static)
112			for (i=0;i<n;++i) mis_marker[i]=-1;
113	artak	1954	/Paso_Pattern_RS(A_p,mis_marker,0.25);/
114			/Paso_Pattern_Aggregiation(A_p,mis_marker,0.5);/
115			Paso_Pattern_coup(A_p,mis_marker,0.05);
116	jfenwick	1851	time2=Paso_timer()-time0;
117			if (Paso_noError()) {
118			#pragma omp parallel for private(i) schedule(static)
119			for (i = 0; i < n; ++i) counter[i]=mis_marker[i];
120			out->n=n;
121			out->n_block=n_block;
122			out->n_F=Paso_Util_cumsum(n,counter);
123			out->mask_F=MEMALLOC(n,index_t);
124			out->rows_in_F=MEMALLOC(out->n_F,index_t);
125			out->inv_A_FF=MEMALLOC(n_blockn_blockout->n_F,double);
126			out->A_FF_pivot=NULL; /* later use for block size>3 */
127			if (! (Paso_checkPtr(out->mask_F) \|\| Paso_checkPtr(out->inv_A_FF) \|\| Paso_checkPtr(out->rows_in_F) ) ) {
128			#pragma omp parallel
129			{
130			/* creates an index for F from mask */
131			#pragma omp for private(i) schedule(static)
132			for (i = 0; i < out->n_F; ++i) out->rows_in_F[i]=-1;
133			#pragma omp for private(i) schedule(static)
134			for (i = 0; i < n; ++i) {
135			if (mis_marker[i]) {
136			out->rows_in_F[counter[i]]=i;
137			out->mask_F[i]=counter[i];
138			} else {
139			out->mask_F[i]=-1;
140			}
141			}
142	artak	1881	/* Compute row-sum for getting rs(A_FF)*/
143			#pragma omp for private(i,iPtr) schedule(static)
144			for (i = 0; i < out->n_F; ++i) {
145	artak	1902	out->inv_A_FF[i]=0;
146	artak	1881	for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) {
147	artak	1902	out->inv_A_FF[i]+=A_p->val[iPtr];
148	artak	1881	}
149			}
150	artak	1954
151			/* for (i=0;i<out->n_F;++i) {
152			fprintf(stderr,"GET INV_A_FF %g %d of %d (%d)\n",out->inv_A_FF[i],i,out->rows_in_F[i],(out->inv_A_FF[i]==0));
153			}
154			fprintf(stderr,"END\n");
155			*/
156	phornby	1917	#pragma omp for private(i, where_p,iPtr,index) schedule(static)
157	jfenwick	1851	for (i = 0; i < out->n_F; i++) {
158			/* find main diagonal */
159			iPtr=A_p->pattern->ptr[out->rows_in_F[i]];
160			index=&(A_p->pattern->index[iPtr]);
161			where_p=(index_t*)bsearch(&out->rows_in_F[i],
162			index,
163			A_p->pattern->ptr[out->rows_in_F[i] + 1]-A_p->pattern->ptr[out->rows_in_F[i]],
164			sizeof(index_t),
165			Paso_comparIndex);
166			if (where_p==NULL) {
167			Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
168			} else {
169			iPtr+=(index_t)(where_p-index);
170			/* get inverse of A_FF block: */
171	artak	1902	S=out->inv_A_FF[i];
172			if (ABS(A_p->val[iPtr])>0.) {
173			if(ABS(S)>0.)
174			out->inv_A_FF[i]=1./S;
175	artak	1954	else
176	artak	1931	{
177	artak	1954	fprintf(stderr,"ROWSUM OF ROW %d->%d is 0\n",i,out->rows_in_F[i]);
178	artak	1931	S=0;
179			for (iPtr=A_p->pattern->ptr[out->rows_in_F[i]];iPtr<A_p->pattern->ptr[out->rows_in_F[i] + 1]; ++iPtr) {
180			if(A_p->val[iPtr]!=0) {
181			fprintf(stderr,"A_p[%d,%d]=%f ",out->rows_in_F[i],A_p->pattern->index[iPtr],A_p->val[iPtr]);
182			S+=A_p->val[iPtr];
183			}
184			}
185			fprintf(stderr,"\n SUMMMMMMM %f\n",S);
186			}
187	artak	1954	} else {
188	artak	1902	Paso_setError(ZERO_DIVISION_ERROR, "Paso_Solver_getAMG: Break-down in AMG decomposition: non-regular main diagonal block.");
189	artak	1890	}
190	artak	1881	}
191	jfenwick	1851	}
192			} /* end parallel region */
193
194			if( Paso_noError()) {
195			/* if there are no nodes in the coarse level there is no more work to do */
196			out->n_C=n-out->n_F;
197	artak	1954	if (level<2) {
198			/if (out->n_C>10) {/
199	jfenwick	1851	out->rows_in_C=MEMALLOC(out->n_C,index_t);
200			out->mask_C=MEMALLOC(n,index_t);
201			if (! (Paso_checkPtr(out->mask_C) \|\| Paso_checkPtr(out->rows_in_C) ) ) {
202			/* creates an index for C from mask */
203			#pragma omp parallel for private(i) schedule(static)
204			for (i = 0; i < n; ++i) counter[i]=! mis_marker[i];
205			Paso_Util_cumsum(n,counter);
206			#pragma omp parallel
207			{
208			#pragma omp for private(i) schedule(static)
209			for (i = 0; i < out->n_C; ++i) out->rows_in_C[i]=-1;
210			#pragma omp for private(i) schedule(static)
211			for (i = 0; i < n; ++i) {
212			if (! mis_marker[i]) {
213			out->rows_in_C[counter[i]]=i;
214			out->mask_C[i]=counter[i];
215			} else {
216			out->mask_C[i]=-1;
217			}
218			}
219			} /* end parallel region */
220			/* get A_CF block: */
221			out->A_CF=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_F,out->rows_in_C,out->mask_F);
222			if (Paso_noError()) {
223			/* get A_FC block: */
224			out->A_FC=Paso_SparseMatrix_getSubmatrix(A_p,out->n_F,out->n_C,out->rows_in_F,out->mask_C);
225	artak	1881	/* get A_CC block: */
226	jfenwick	1851	if (Paso_noError()) {
227			schur=Paso_SparseMatrix_getSubmatrix(A_p,out->n_C,out->n_C,out->rows_in_C,out->mask_C);
228	artak	1881
229			/find the pattern of the schur complement with fill in/
230			schur_withFillIn=Paso_SparseMatrix_alloc(A_p->type,Paso_Pattern_binop(PATTERN_FORMAT_DEFAULT, schur->pattern, Paso_Pattern_multiply(PATTERN_FORMAT_DEFAULT,out->A_CF->pattern,out->A_FC->pattern)),1,1);
231	artak	1902
232	artak	1881	/* copy values over*/
233			#pragma omp for private(i,iPtr,iPtr_s,j,j0) schedule(static)
234			for (i = 0; i < schur_withFillIn->numRows; ++i) {
235			for (iPtr=schur_withFillIn->pattern->ptr[i];iPtr<schur_withFillIn->pattern->ptr[i + 1]; ++iPtr) {
236			j=schur_withFillIn->pattern->index[iPtr];
237			schur_withFillIn->val[iPtr]=0.;
238			for (iPtr_s=schur->pattern->ptr[i];iPtr_s<schur->pattern->ptr[i + 1]; ++iPtr_s){
239			j0=schur->pattern->index[iPtr_s];
240			if (j==j0) {
241			schur_withFillIn->val[iPtr]=schur->val[iPtr_s];
242			break;
243			}
244			}
245			}
246			}
247	artak	1902	time0=Paso_timer()-time0;
248	artak	1931
249	jfenwick	1851	if (Paso_noError()) {
250			time1=Paso_timer();
251			/* update A_CC block to get Schur complement and then apply AMG to it */
252	artak	1881	Paso_Solver_updateIncompleteSchurComplement(schur_withFillIn,out->A_CF,out->inv_A_FF,out->A_FF_pivot,out->A_FC);
253	jfenwick	1851	time1=Paso_timer()-time1;
254	artak	1939	out->AMG_of_Schur=Paso_Solver_getAMG(schur_withFillIn,verbose,level+1);
255	artak	1881
256	artak	1931	/Paso_SparseMatrix_free(schur);/
257	artak	1890	/* Paso_SparseMatrix_free(schur_withFillIn);*/
258	jfenwick	1851	}
259			/* allocate work arrays for AMG application */
260			if (Paso_noError()) {
261			out->x_F=MEMALLOC(n_block*out->n_F,double);
262			out->b_F=MEMALLOC(n_block*out->n_F,double);
263			out->x_C=MEMALLOC(n_block*out->n_C,double);
264			out->b_C=MEMALLOC(n_block*out->n_C,double);
265	artak	1890
266	jfenwick	1851	if (! (Paso_checkPtr(out->x_F) \|\| Paso_checkPtr(out->b_F) \|\| Paso_checkPtr(out->x_C) \|\| Paso_checkPtr(out->b_C) ) ) {
267			#pragma omp parallel
268			{
269			#pragma omp for private(i,k) schedule(static)
270			for (i = 0; i < out->n_F; ++i) {
271			for (k=0; k<n_block;++k) {
272			out->x_F[i*n_block+k]=0.;
273			out->b_F[i*n_block+k]=0.;
274			}
275			}
276			#pragma omp for private(i,k) schedule(static)
277			for (i = 0; i < out->n_C; ++i) {
278			for (k=0; k<n_block;++k) {
279			out->x_C[i*n_block+k]=0.;
280			out->b_C[i*n_block+k]=0.;
281			}
282			}
283			} /* end parallel region */
284			}
285			}
286			}
287			}
288			}
289			}
290			}
291			}
292			}
293			}
294			TMPMEMFREE(mis_marker);
295			TMPMEMFREE(counter);
296			if (Paso_noError()) {
297			if (verbose) {
298			printf("AMG: %d unknowns eliminated. %d left.\n",out->n_F,n-out->n_F);
299	artak	1954	if (out->n_C>10) {
300	jfenwick	1851	printf("timing: AMG: MIS/reordering/elemination : %e/%e/%e\n",time2,time0,time1);
301			} else {
302			printf("timing: AMG: MIS: %e\n",time2);
303			}
304			}
305			return out;
306			} else {
307			Paso_Solver_AMG_free(out);
308			return NULL;
309			}
310			}
311
312			/**************************************************************/
313
314			/* apply AMG precondition b-> x
315
316			in fact it solves
317
318			[ I 0 ] [ A_FF 0 ] [ I invA_FF*A_FF ] [ x_F ] = [b_F]
319			[ A_CF*invA_FF I ] [ 0 S ] [ 0 I ] [ x_C ] = [b_C]
320
321			in the form
322
323			b->[b_F,b_C]
324			x_F=invA_FF*b_F
325			b_C=b_C-A_CF*x_F
326			x_C=AMG(b_C)
327			b_F=b_F-A_FC*x_C
328			x_F=invA_FF*b_F
329			x<-[x_F,x_C]
330
331			should be called within a parallel region
332			barrier synconization should be performed to make sure that the input vector available
333
334			*/
335
336			void Paso_Solver_solveAMG(Paso_Solver_AMG * amg, double * x, double * b) {
337	artak	1954	dim_t i,k,oldsweeps,j;
338	artak	1890	double *r=MEMALLOC(amg->n,double);
339	artak	1939	/Paso_Solver_GS GS=NULL;*/
340	artak	1954	double *bold=MEMALLOC(amg->n,double);
341			double *bnew=MEMALLOC(amg->n,double);
342			double *x0=MEMALLOC(amg->n,double);
343			index_t iPtr;
344	artak	1939
345	artak	1954	if (amg->level==2) {
346			/if (amg->n_C<=10) {/
347	artak	1939	Paso_UMFPACK1(amg->A,x,b,0);
348	jfenwick	1851	} else {
349	artak	1954
350	artak	1902	/* presmoothing */
351	artak	1939	Paso_Solver_solveGS(amg->GS,x,b);
352			oldsweeps=amg->GS->sweeps;
353			if (amg->GS->sweeps>1) {
354	artak	1931
355			#pragma omp parallel for private(i) schedule(static)
356	artak	1954	for (i=0;i<amg->GS->n;++i) bold[i]=b[i];
357	artak	1931
358	artak	1939	while(amg->GS->sweeps>1) {
359	artak	1931	#pragma omp parallel for private(i) schedule(static)
360	artak	1954	for (i=0;i<amg->GS->n;++i) bnew[i]=bold[i]+b[i];
361	artak	1931	/* Compute the residual b=b-Ax*/
362			Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(DBLE(-1), amg->A, x, DBLE(1), bnew);
363			/* Go round again*/
364	artak	1939	Paso_Solver_solveGS(amg->GS,x,bnew);
365	artak	1931	#pragma omp parallel for private(i) schedule(static)
366	artak	1954	for (i=0;i<amg->GS->n;++i) bold[i]=bnew[i];
367	artak	1939	amg->GS->sweeps=amg->GS->sweeps-1;
368	artak	1931	}
369			}
370	artak	1939	amg->GS->sweeps=oldsweeps;
371	artak	1931	/* end of presmoothing */
372	artak	1939
373	artak	1890	#pragma omp parallel for private(i) schedule(static)
374			for (i=0;i<amg->n;++i) r[i]=b[i];
375
376			/r=b-Ax/
377			Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
378	artak	1954
379	jfenwick	1851	/* b->[b_F,b_C] */
380	artak	1954	#pragma omp parallel for private(i) schedule(static)
381			for (i=0;i<amg->n_F;++i) amg->b_F[i]=r[amg->rows_in_F[i]];
382
383			#pragma omp parallel for private(i) schedule(static)
384			for (i=0;i<amg->n_C;++i) amg->b_C[i]=r[amg->rows_in_C[i]];
385	artak	1890
386	jfenwick	1851	/* x_F=invA_FFb_F /
387	artak	1954	Paso_Solver_applyBlockDiagonalMatrix(1,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F);
388	artak	1890
389	jfenwick	1851	/* b_C=b_C-A_CFx_F /
390			Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_CF,amg->x_F,1.,amg->b_C);
391	artak	1890
392	jfenwick	1851	/* x_C=AMG(b_C) */
393			Paso_Solver_solveAMG(amg->AMG_of_Schur,amg->x_C,amg->b_C);
394	artak	1890
395	jfenwick	1851	/* b_F=b_F-A_FCx_C /
396			Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A_FC,amg->x_C,1.,amg->b_F);
397			/* x_F=invA_FFb_F /
398	artak	1954	Paso_Solver_applyBlockDiagonalMatrix(1,amg->n_F,amg->inv_A_FF,amg->A_FF_pivot,amg->x_F,amg->b_F);
399	jfenwick	1851	/* x<-[x_F,x_C] */
400	artak	1890
401	artak	1954	#pragma omp parallel for private(i) schedule(static)
402			for (i=0;i<amg->n;++i) {
403			if (amg->mask_C[i]>-1) {
404			x[i]+=amg->x_C[amg->mask_C[i]];
405			} else {
406			x[i]+=amg->x_F[amg->mask_F[i]];
407			}
408	jfenwick	1851	}
409	artak	1954
410	artak	1931
411			/postsmoothing/
412	artak	1954	#pragma omp parallel for private(i) schedule(static)
413			for (i=0;i<amg->n;++i) r[i]=b[i];
414
415			/for (i=0;i<amg->n;++i) fprintf(stderr,"%lf %lf %d %d %d BBB\n",r[i],x[i],amg->n,amg->A->numCols, amg->level);/
416
417			/*
418			for (i = 0; i < amg->A->numRows; ++i) {
419			for (iPtr=amg->A->pattern->ptr[i];iPtr<amg->A->pattern->ptr[i + 1]; ++iPtr) {
420			j=amg->A->pattern->index[iPtr];
421			fprintf(stderr,"A[%d,%d]=%lf ",i,j,amg->A->val[iPtr]);
422			}
423			fprintf(stderr,"\n");
424			}*/
425
426			/r=b-Ax/
427			Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(-1.,amg->A,x,1.,r);
428			Paso_Solver_solveGS(amg->GS,x0,r);
429
430			#pragma omp parallel for private(i) schedule(static)
431			for (i=0;i<amg->n;++i) {
432			/fprintf(stderr,"%lf %lf %lf\n",x[i],x0[i],x[i]+x0[i]);/
433			x[i]+=x0[i];
434			}
435
436	artak	1931	/end of postsmoothing/
437
438	artak	1954	/Paso_Solver_GS_free(amg->GS);/
439	jfenwick	1851	}
440	artak	1890	MEMFREE(r);
441	artak	1931	MEMFREE(bold);
442	artak	1954	MEMFREE(bnew);
443			MEMFREE(x0);
444	artak	1931	return;
445	jfenwick	1851	}
446
447			/*
448			* $Log$
449			*
450			*/