paso/src/Pattern_coupling.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: Pattern: Paso_Pattern_coupling 

   searches for a maximal independent set MIS in the matrix pattern 
   vertices in the maximal independent set are marked in mis_marker
   nodes to be considered are marked by -1 on the input in mis_marker

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

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

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

#include "PasoUtil.h"
#include "Pattern_coupling.h"
#include <limits.h>


/***************************************************************/
 
#define IS_AVAILABLE -1
#define IS_IN_SET -3
#define IS_REMOVED -4

void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {

  dim_t i,j;
  /*double threshold=0.05;*/
  double sum;
  index_t iptr,*index,*where_p,*diagptr;
  bool_t passed=FALSE;
  dim_t n=A->numRows;
  diagptr=MEMALLOC(n,index_t);

  if (A->pattern->type & PATTERN_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
    return;
  }
   
   #pragma omp parallel for private(i) schedule(static)
   for (i=0;i<n;++i)
        if(mis_marker[i]==IS_AVAILABLE)
                    mis_marker[i]=IS_IN_SET;

    #pragma omp parallel for private(i,index,where_p) schedule(static) 
    for (i=0;i<n;++i) {
         diagptr[i]=A->pattern->ptr[i];
         index=&(A->pattern->index[diagptr[i]]);
         where_p=(index_t*)bsearch(&i,
                                index,
                                A->pattern->ptr[i + 1]-A->pattern->ptr[i],
                                sizeof(index_t),
                                Paso_comparIndex);
        if (where_p==NULL) {
            Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
        } else {
                diagptr[i]+=(index_t)(where_p-index);
        }
    }
    

    /*This loop cannot be parallelized, as order matters here.*/ 
    for (i=0;i<n;++i) {
      if (mis_marker[i]==IS_IN_SET) {
        for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
             j=A->pattern->index[iptr];
             if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
                mis_marker[j]=IS_REMOVED;
             }
        }
      }
    }
    
    
      /*This loop cannot be parallelized, as order matters here.*/ 
    for (i=0;i<n;i++) {
        if (mis_marker[i]==IS_REMOVED) {
           for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
              j=A->pattern->index[iptr];
              if (mis_marker[j]==IS_IN_SET) {
                if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
                    passed=TRUE;
                }
                else {
                    passed=FALSE;
                    break;
                }
              } 
            }
            if (passed) {
               mis_marker[i]=IS_IN_SET;
               passed=FALSE;
            }
        }
    }

    /* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
    /* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
    #pragma omp parallel for private(i,iptr,j,sum) schedule(static)
    for (i=0;i<n;i++) {
        if (mis_marker[i]==IS_REMOVED) {
           sum=0;
           for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
             j=A->pattern->index[iptr];
             if (mis_marker[j]==IS_REMOVED)
                sum+=A->val[iptr];
           }
           if(ABS(sum)<1.e-25)
             mis_marker[i]=IS_IN_SET;
        }
    }


     /* swap to TRUE/FALSE in mis_marker */
     #pragma omp parallel for private(i) schedule(static)
     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
     
     MEMFREE(diagptr);
}

/*
 * Ruge-Stueben strength of connection mask.
 *
 */
void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
{
  dim_t i,n,j;
  index_t iptr;
  double threshold,min_offdiagonal;
  
  Paso_Pattern *out=NULL;
  
  Paso_IndexList* index_list=NULL;

 index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
   if (! Paso_checkPtr(index_list)) {
        #pragma omp parallel for private(i) schedule(static)
        for(i=0;i<A->pattern->numOutput;++i) {
             index_list[i].extension=NULL;
             index_list[i].n=0;
        }
    }
  
  
  n=A->numRows;
  if (A->pattern->type & PATTERN_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
    return;
  }

    #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
    for (i=0;i<n;++i) {
      if(mis_marker[i]==IS_AVAILABLE) {
        min_offdiagonal = DBL_MAX;
        for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
            if(A->pattern->index[iptr] != i){
                min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
            }
        }

        threshold = theta*min_offdiagonal;
        for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
            j=A->pattern->index[iptr];
            if(A->val[iptr]<=threshold) {
               if(j!=i) {
                Paso_IndexList_insertIndex(&(index_list[i]),j);
                }
            }
        }
       }
      }
    
    out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
    
     /* clean up */
   if (index_list!=NULL) {
        #pragma omp parallel for private(i) schedule(static)
        for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
     }
    TMPMEMFREE(index_list);

    Paso_Pattern_mis(out,mis_marker);
}

void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
{
  dim_t i,j,n;
  index_t iptr;
  double diag,eps_Aii,val;
  double* diags;


  Paso_Pattern *out=NULL;
  Paso_IndexList* index_list=NULL;

  n=A->numRows;  
  diags=MEMALLOC(n,double);

  index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
   if (! Paso_checkPtr(index_list)) {
        #pragma omp parallel for private(i) schedule(static)
        for(i=0;i<A->pattern->numOutput;++i) {
             index_list[i].extension=NULL;
             index_list[i].n=0;
        }
    }
    
  if (A->pattern->type & PATTERN_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
    return;
  }


    #pragma omp parallel for private(i,iptr,diag) schedule(static)
      for (i=0;i<n;++i) {
        diag = 0;
        for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
            if(A->pattern->index[iptr] != i){
                diag+=A->val[iptr];
            }
        }
        diags[i]=ABS(diag);
      }


    #pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
     for (i=0;i<n;++i) {
       if (mis_marker[i]==IS_AVAILABLE) {
        eps_Aii = theta*theta*diags[i];
        val=0.;
        for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
            j=A->pattern->index[iptr];
            val=A->val[iptr];
            if(j!= i) {
              if(val*val>=eps_Aii * diags[j]) {
               Paso_IndexList_insertIndex(&(index_list[i]),j);
              }
            }
        }
       }
     }

    out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
    
     /* clean up */
    if (index_list!=NULL) {
        #pragma omp parallel for private(i) schedule(static)
        for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
     }

    TMPMEMFREE(index_list);
    MEMFREE(diags);
    
    Paso_Pattern_mis(out,mis_marker);


}


#undef IS_AVAILABLE 
#undef IS_IN_SET 
#undef IS_REMOVED
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: Pattern: Paso_Pattern_coupling
18
19			searches for a maximal independent set MIS in the matrix pattern
20			vertices in the maximal independent set are marked in mis_marker
21			nodes to be considered are marked by -1 on the input in mis_marker
22
23			*/
24			/**********************************************************************/
25
26			/* Copyrights by ACcESS Australia 2003,2004,2005 */
27			/* Author: artak@uq.edu.au */
28
29			/**************************************************************/
30
31	phornby	1914	#include "PasoUtil.h"
32	phornby	1913	#include "Pattern_coupling.h"
33	artak	2381	#include <limits.h>
34	jfenwick	1851
35
36			/***************************************************************/
37
38			#define IS_AVAILABLE -1
39	artak	1975	#define IS_IN_SET -3
40			#define IS_REMOVED -4
41	jfenwick	1851
42	artak	1931	void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43	jfenwick	1851
44	artak	1881	dim_t i,j;
45	artak	1931	/double threshold=0.05;/
46	artak	2122	double sum;
47	artak	2107	index_t iptr,index,where_p,*diagptr;
48			bool_t passed=FALSE;
49	artak	1954	dim_t n=A->numRows;
50	artak	2107	diagptr=MEMALLOC(n,index_t);
51	artak	2122
52	jfenwick	1851	if (A->pattern->type & PATTERN_FORMAT_SYM) {
53	artak	2159	Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
54	jfenwick	1851	return;
55			}
56
57	artak	2107	#pragma omp parallel for private(i) schedule(static)
58			for (i=0;i<n;++i)
59			if(mis_marker[i]==IS_AVAILABLE)
60			mis_marker[i]=IS_IN_SET;
61	jfenwick	1851
62	artak	2107	#pragma omp parallel for private(i,index,where_p) schedule(static)
63			for (i=0;i<n;++i) {
64			diagptr[i]=A->pattern->ptr[i];
65			index=&(A->pattern->index[diagptr[i]]);
66			where_p=(index_t*)bsearch(&i,
67			index,
68			A->pattern->ptr[i + 1]-A->pattern->ptr[i],
69			sizeof(index_t),
70			Paso_comparIndex);
71			if (where_p==NULL) {
72	artak	2159	Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
73	artak	2107	} else {
74			diagptr[i]+=(index_t)(where_p-index);
75			}
76			}
77
78
79
80			/This loop cannot be parallelized, as order matters here./
81			for (i=0;i<n;++i) {
82			if (mis_marker[i]==IS_IN_SET) {
83			for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
84			j=A->pattern->index[iptr];
85			if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
86			mis_marker[j]=IS_REMOVED;
87			}
88			}
89			}
90			}
91
92
93
94			/This loop cannot be parallelized, as order matters here./
95			for (i=0;i<n;i++) {
96			if (mis_marker[i]==IS_REMOVED) {
97			for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
98			j=A->pattern->index[iptr];
99			if (mis_marker[j]==IS_IN_SET) {
100			if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
101			passed=TRUE;
102	artak	1902	}
103	artak	2107	else {
104			passed=FALSE;
105			break;
106	artak	1902	}
107	artak	2307	}
108	artak	1902	}
109	artak	2107	if (passed) {
110			mis_marker[i]=IS_IN_SET;
111			passed=FALSE;
112			}
113	artak	1902	}
114	artak	2107	}
115
116	artak	2381	/* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
117	artak	2122	/* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
118	caltinay	2309	#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
119	artak	2122	for (i=0;i<n;i++) {
120			if (mis_marker[i]==IS_REMOVED) {
121			sum=0;
122			for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
123			j=A->pattern->index[iptr];
124			if (mis_marker[j]==IS_REMOVED)
125			sum+=A->val[iptr];
126			}
127	artak	2381	if(ABS(sum)<1.e-25)
128	artak	2122	mis_marker[i]=IS_IN_SET;
129			}
130			}
131
132
133	jfenwick	1851	/* swap to TRUE/FALSE in mis_marker */
134			#pragma omp parallel for private(i) schedule(static)
135	artak	2107	for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
136
137			MEMFREE(diagptr);
138	jfenwick	1851	}
139	artak	1890
140	artak	2381	/*
141			* Ruge-Stueben strength of connection mask.
142			*
143			*/
144			void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
145			{
146			dim_t i,n,j;
147			index_t iptr;
148			double threshold,min_offdiagonal;
149
150			Paso_Pattern *out=NULL;
151
152			Paso_IndexList* index_list=NULL;
153	artak	1890
154	artak	2381	index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
155			if (! Paso_checkPtr(index_list)) {
156			#pragma omp parallel for private(i) schedule(static)
157			for(i=0;i<A->pattern->numOutput;++i) {
158			index_list[i].extension=NULL;
159			index_list[i].n=0;
160			}
161			}
162
163
164			n=A->numRows;
165			if (A->pattern->type & PATTERN_FORMAT_SYM) {
166			Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
167			return;
168			}
169
170			#pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
171			for (i=0;i<n;++i) {
172			if(mis_marker[i]==IS_AVAILABLE) {
173			min_offdiagonal = DBL_MAX;
174			for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
175			if(A->pattern->index[iptr] != i){
176			min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
177			}
178			}
179
180			threshold = theta*min_offdiagonal;
181			for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
182			j=A->pattern->index[iptr];
183			if(A->val[iptr]<=threshold) {
184			if(j!=i) {
185			Paso_IndexList_insertIndex(&(index_list[i]),j);
186			}
187			}
188			}
189			}
190			}
191
192			out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
193
194			/* clean up */
195			if (index_list!=NULL) {
196			#pragma omp parallel for private(i) schedule(static)
197			for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
198			}
199			TMPMEMFREE(index_list);
200
201			Paso_Pattern_mis(out,mis_marker);
202			}
203
204			void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
205			{
206			dim_t i,j,n;
207			index_t iptr;
208			double diag,eps_Aii,val;
209			double* diags;
210
211
212			Paso_Pattern *out=NULL;
213			Paso_IndexList* index_list=NULL;
214
215			n=A->numRows;
216			diags=MEMALLOC(n,double);
217
218			index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
219			if (! Paso_checkPtr(index_list)) {
220			#pragma omp parallel for private(i) schedule(static)
221			for(i=0;i<A->pattern->numOutput;++i) {
222			index_list[i].extension=NULL;
223			index_list[i].n=0;
224			}
225			}
226
227			if (A->pattern->type & PATTERN_FORMAT_SYM) {
228			Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
229			return;
230			}
231
232
233			#pragma omp parallel for private(i,iptr,diag) schedule(static)
234			for (i=0;i<n;++i) {
235			diag = 0;
236			for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
237			if(A->pattern->index[iptr] != i){
238			diag+=A->val[iptr];
239			}
240			}
241			diags[i]=ABS(diag);
242			}
243
244
245			#pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
246			for (i=0;i<n;++i) {
247			if (mis_marker[i]==IS_AVAILABLE) {
248			eps_Aii = thetathetadiags[i];
249			val=0.;
250			for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
251			j=A->pattern->index[iptr];
252			val=A->val[iptr];
253			if(j!= i) {
254			if(valval>=eps_Aii diags[j]) {
255			Paso_IndexList_insertIndex(&(index_list[i]),j);
256			}
257			}
258			}
259			}
260			}
261
262			out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
263
264			/* clean up */
265			if (index_list!=NULL) {
266			#pragma omp parallel for private(i) schedule(static)
267			for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
268			}
269
270			TMPMEMFREE(index_list);
271			MEMFREE(diags);
272
273			Paso_Pattern_mis(out,mis_marker);
274
275
276			}
277
278
279
280
281
282	jfenwick	1851	#undef IS_AVAILABLE
283	artak	1975	#undef IS_IN_SET
284			#undef IS_REMOVED