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
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	#include "PasoUtil.h"
32	#include "Pattern_coupling.h"
33	#include <limits.h>
34
35
36	/***************************************************************/
37
38	#define IS_AVAILABLE -1
39	#define IS_IN_SET -3
40	#define IS_REMOVED -4
41
42	void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43
44	dim_t i,j;
45	/double threshold=0.05;/
46	double sum;
47	index_t iptr,index,where_p,*diagptr;
48	bool_t passed=FALSE;
49	dim_t n=A->numRows;
50	diagptr=MEMALLOC(n,index_t);
51
52	if (A->pattern->type & PATTERN_FORMAT_SYM) {
53	Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
54	return;
55	}
56
57	#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
62	#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	Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
73	} 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	}
103	else {
104	passed=FALSE;
105	break;
106	}
107	}
108	}
109	if (passed) {
110	mis_marker[i]=IS_IN_SET;
111	passed=FALSE;
112	}
113	}
114	}
115
116	/* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
117	/* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
118	#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
119	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	if(ABS(sum)<1.e-25)
128	mis_marker[i]=IS_IN_SET;
129	}
130	}
131
132
133	/* swap to TRUE/FALSE in mis_marker */
134	#pragma omp parallel for private(i) schedule(static)
135	for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
136
137	MEMFREE(diagptr);
138	}
139
140	/*
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
154	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	#undef IS_AVAILABLE
283	#undef IS_IN_SET
284	#undef IS_REMOVED