paso/src/Pattern_coupling.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: 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   /* Week connection */
#define IS_REMOVED -4  /* strong */

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

  dim_t i,j;
  /*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_REMOVED;

    #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[A->pattern->ptr[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_REMOVED) {
        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_IN_SET;
             }
        }
      }
    }
    
    
      /*This loop cannot be parallelized, as order matters here.*/ 
    for (i=0;i<n;i++) {
        if (mis_marker[i]==IS_IN_SET) {
           passed=TRUE;
           for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
              j=A->pattern->index[iptr];
              if (mis_marker[j]==IS_REMOVED) {
                if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
                    passed=TRUE;
                }
                else {
                    passed=FALSE;
                    break;
                }
              } 
           }
           if (passed) mis_marker[i]=IS_REMOVED;
        }
    }
    /* 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,max_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,max_offdiagonal,threshold,j) schedule(static)*/
    for (i=0;i<n;++i) {
      if(mis_marker[i]==IS_AVAILABLE) {
        max_offdiagonal = DBL_MIN;
        for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
            if(A->pattern->index[iptr] != i){
                max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
            }
        }
        
        threshold = theta*max_offdiagonal;
        for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
            j=A->pattern->index[iptr];
            if((-A->val[iptr])>=threshold) {
                Paso_IndexList_insertIndex(&(index_list[i]),j);
                Paso_IndexList_insertIndex(&(index_list[j]),i);
            }
        }
       }
      }
    
   
    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);*/
    Paso_Pattern_greedy(out,mis_marker);
    Paso_Pattern_free(out);
}

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_Aggregiation: 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];
        for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
            j=A->pattern->index[iptr];
            val=A->val[iptr];
              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);*/
    Paso_Pattern_greedy(out,mis_marker);
    Paso_Pattern_free(out);

}

/* Greedy algorithm */
void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {

  dim_t i,j;
  /*double sum;*/
  index_t iptr;
  bool_t passed=FALSE;
  dim_t n=pattern->numOutput;

  if (pattern->type & PATTERN_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
    return;
  }
   
   /* We do not need this loop if we set IS_REMOVED=IS_AVAILABLE. */
   #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;


    for (i=0;i<n;++i) {
      if (mis_marker[i]==IS_IN_SET) {
        for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
             j=pattern->index[iptr];
             mis_marker[j]=IS_REMOVED;
        }
      }
    }
    
    
    for (i=0;i<n;i++) {
        if (mis_marker[i]==IS_REMOVED) {
           passed=TRUE;
           for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
              j=pattern->index[iptr];
                if (mis_marker[j]==IS_REMOVED) {
                    passed=TRUE;
                }
                else {
                    passed=FALSE;
                    break;
                }
              } 
           if (passed) 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_REMOVED);
     
}


void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {

  dim_t i,j;
  /*double sum;*/
  index_t iptr;
  index_t num_colors;
  index_t* colorOf;
  register index_t color;
  bool_t passed=FALSE;
  dim_t n=pattern->numOutput;

  
  colorOf=MEMALLOC(n,index_t);

  if (pattern->type & PATTERN_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
    return;
  }
   
   Paso_Pattern_color(pattern,&num_colors,colorOf);
   
   /* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
   #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 barrier
   for (color=0;color<num_colors;++color) {
    #pragma omp parallel for schedule(static) private(i,iptr,j)
    for (i=0;i<n;++i) {
     if (colorOf[i]==color) {  
      if (mis_marker[i]==IS_IN_SET) {
        for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
             j=pattern->index[iptr];
             if (colorOf[j]<color)
              mis_marker[j]=IS_REMOVED;
        }
      }
     }
    }
   }
    
    
   #pragma omp barrier
   for (color=0;color<num_colors;++color) {
   #pragma omp parallel for schedule(static) private(i,iptr,j) 
    for (i=0;i<n;i++) {
      if (colorOf[i]==color) {  
        if (mis_marker[i]==IS_REMOVED) {
           passed=TRUE;
           for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
              j=pattern->index[iptr];
               if (colorOf[j]<color && passed) {
                if (mis_marker[j]==IS_REMOVED) {
                    passed=TRUE;
                }
                else {
                    passed=FALSE;
                    /*break;*/
                }
              }
           }
           if (passed) 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(colorOf); 
}

#undef IS_AVAILABLE 
#undef IS_IN_SET 
#undef IS_REMOVED


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: 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 /* Week connection */
40	#define IS_REMOVED -4 /* strong */
41
42	void Paso_Pattern_YS(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
43
44	dim_t i,j;
45	/double sum;/
46	index_t iptr,index,where_p,*diagptr;
47	bool_t passed=FALSE;
48	dim_t n=A->numRows;
49	diagptr=MEMALLOC(n,index_t);
50
51	if (A->pattern->type & PATTERN_FORMAT_SYM) {
52	Paso_setError(TYPE_ERROR,"Paso_Pattern_coup: symmetric matrix pattern is not supported yet");
53	return;
54	}
55
56	#pragma omp parallel for private(i) schedule(static)
57	for (i=0;i<n;++i)
58	if(mis_marker[i]==IS_AVAILABLE)
59	mis_marker[i]=IS_REMOVED;
60
61	#pragma omp parallel for private(i,index,where_p) schedule(static)
62	for (i=0;i<n;++i) {
63	diagptr[i]=A->pattern->ptr[i];
64	index=&(A->pattern->index[A->pattern->ptr[i]]);
65	where_p=(index_t*)bsearch(&i,
66	index,
67	A->pattern->ptr[i + 1]-A->pattern->ptr[i],
68	sizeof(index_t),
69	Paso_comparIndex);
70	if (where_p==NULL) {
71	Paso_setError(VALUE_ERROR, "Paso_Pattern_coup: main diagonal element missing.");
72	} else {
73	diagptr[i]+=(index_t)(where_p-index);
74	}
75	}
76
77
78
79	/This loop cannot be parallelized, as order matters here./
80	for (i=0;i<n;++i) {
81	if (mis_marker[i]==IS_REMOVED) {
82	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
83	j=A->pattern->index[iptr];
84	if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
85	mis_marker[j]=IS_IN_SET;
86	}
87	}
88	}
89	}
90
91
92
93	/This loop cannot be parallelized, as order matters here./
94	for (i=0;i<n;i++) {
95	if (mis_marker[i]==IS_IN_SET) {
96	passed=TRUE;
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_REMOVED) {
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) mis_marker[i]=IS_REMOVED;
110	}
111	}
112	/* This check is to make sure we dont get some nusty rows which were not removed durring coarsening process.*/
113	/* TODO: we have to mechanism that this does not happend at all, and get rid of this 'If'. */
114	/*#pragma omp parallel for private(i,iptr,j,sum) schedule(static)
115	for (i=0;i<n;i++) {
116	if (mis_marker[i]==IS_REMOVED) {
117	sum=0;
118	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
119	j=A->pattern->index[iptr];
120	if (mis_marker[j]==IS_REMOVED)
121	sum+=A->val[iptr];
122	}
123	if(ABS(sum)<1.e-25)
124	mis_marker[i]=IS_IN_SET;
125	}
126	}
127	*/
128
129	/* swap to TRUE/FALSE in mis_marker */
130	#pragma omp parallel for private(i) schedule(static)
131	for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
132
133	MEMFREE(diagptr);
134	}
135
136	/*
137	* Ruge-Stueben strength of connection mask.
138	*
139	*/
140	void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
141	{
142	dim_t i,n,j;
143	index_t iptr;
144	double threshold,max_offdiagonal;
145
146	Paso_Pattern *out=NULL;
147
148	Paso_IndexList* index_list=NULL;
149
150	index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
151	if (! Paso_checkPtr(index_list)) {
152	#pragma omp parallel for private(i) schedule(static)
153	for(i=0;i<A->pattern->numOutput;++i) {
154	index_list[i].extension=NULL;
155	index_list[i].n=0;
156	}
157	}
158
159
160	n=A->numRows;
161	if (A->pattern->type & PATTERN_FORMAT_SYM) {
162	Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
163	return;
164	}
165	/#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold,j) schedule(static)/
166	for (i=0;i<n;++i) {
167	if(mis_marker[i]==IS_AVAILABLE) {
168	max_offdiagonal = DBL_MIN;
169	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
170	if(A->pattern->index[iptr] != i){
171	max_offdiagonal = MAX(max_offdiagonal,-A->val[iptr]);
172	}
173	}
174
175	threshold = theta*max_offdiagonal;
176	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
177	j=A->pattern->index[iptr];
178	if((-A->val[iptr])>=threshold) {
179	Paso_IndexList_insertIndex(&(index_list[i]),j);
180	Paso_IndexList_insertIndex(&(index_list[j]),i);
181	}
182	}
183	}
184	}
185
186
187	out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
188
189	/* clean up */
190	if (index_list!=NULL) {
191	#pragma omp parallel for private(i) schedule(static)
192	for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
193	}
194	TMPMEMFREE(index_list);
195
196	/Paso_Pattern_mis(out,mis_marker);/
197	Paso_Pattern_greedy(out,mis_marker);
198	Paso_Pattern_free(out);
199	}
200
201	void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
202	{
203	dim_t i,j,n;
204	index_t iptr;
205	double diag,eps_Aii,val;
206	double* diags;
207
208
209	Paso_Pattern *out=NULL;
210	Paso_IndexList* index_list=NULL;
211
212	n=A->numRows;
213	diags=MEMALLOC(n,double);
214
215	index_list=TMPMEMALLOC(A->pattern->numOutput,Paso_IndexList);
216	if (! Paso_checkPtr(index_list)) {
217	#pragma omp parallel for private(i) schedule(static)
218	for(i=0;i<A->pattern->numOutput;++i) {
219	index_list[i].extension=NULL;
220	index_list[i].n=0;
221	}
222	}
223
224	if (A->pattern->type & PATTERN_FORMAT_SYM) {
225	Paso_setError(TYPE_ERROR,"Paso_Pattern_Aggregiation: symmetric matrix pattern is not supported yet");
226	return;
227	}
228
229
230	#pragma omp parallel for private(i,iptr,diag) schedule(static)
231	for (i=0;i<n;++i) {
232	diag = 0;
233	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
234	if(A->pattern->index[iptr] == i){
235	diag+=A->val[iptr];
236	}
237	}
238	diags[i]=ABS(diag);
239	}
240
241
242	#pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
243	for (i=0;i<n;++i) {
244	if (mis_marker[i]==IS_AVAILABLE) {
245	eps_Aii = thetathetadiags[i];
246	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
247	j=A->pattern->index[iptr];
248	val=A->val[iptr];
249	if((valval)>=(eps_Aiidiags[j])) {
250	Paso_IndexList_insertIndex(&(index_list[i]),j);
251	}
252	}
253	}
254	}
255
256	out=Paso_IndexList_createPattern(0, A->pattern->numOutput,index_list,0,A->pattern->numInput,0);
257
258	/* clean up */
259	if (index_list!=NULL) {
260	#pragma omp parallel for private(i) schedule(static)
261	for(i=0;i<A->pattern->numOutput;++i) Paso_IndexList_free(index_list[i].extension);
262	}
263
264	TMPMEMFREE(index_list);
265	MEMFREE(diags);
266
267
268	/Paso_Pattern_mis(out,mis_marker);/
269	Paso_Pattern_greedy(out,mis_marker);
270	Paso_Pattern_free(out);
271
272	}
273
274	/* Greedy algorithm */
275	void Paso_Pattern_greedy(Paso_Pattern* pattern, index_t* mis_marker) {
276
277	dim_t i,j;
278	/double sum;/
279	index_t iptr;
280	bool_t passed=FALSE;
281	dim_t n=pattern->numOutput;
282
283	if (pattern->type & PATTERN_FORMAT_SYM) {
284	Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
285	return;
286	}
287
288	/* We do not need this loop if we set IS_REMOVED=IS_AVAILABLE. */
289	#pragma omp parallel for private(i) schedule(static)
290	for (i=0;i<n;++i)
291	if(mis_marker[i]==IS_AVAILABLE)
292	mis_marker[i]=IS_IN_SET;
293
294
295	for (i=0;i<n;++i) {
296	if (mis_marker[i]==IS_IN_SET) {
297	for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
298	j=pattern->index[iptr];
299	mis_marker[j]=IS_REMOVED;
300	}
301	}
302	}
303
304
305
306	for (i=0;i<n;i++) {
307	if (mis_marker[i]==IS_REMOVED) {
308	passed=TRUE;
309	for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
310	j=pattern->index[iptr];
311	if (mis_marker[j]==IS_REMOVED) {
312	passed=TRUE;
313	}
314	else {
315	passed=FALSE;
316	break;
317	}
318	}
319	if (passed) mis_marker[i]=IS_IN_SET;
320	}
321	}
322
323	/* swap to TRUE/FALSE in mis_marker */
324	#pragma omp parallel for private(i) schedule(static)
325	for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_REMOVED);
326
327	}
328
329
330	void Paso_Pattern_greedy_color(Paso_Pattern* pattern, index_t* mis_marker) {
331
332	dim_t i,j;
333	/double sum;/
334	index_t iptr;
335	index_t num_colors;
336	index_t* colorOf;
337	register index_t color;
338	bool_t passed=FALSE;
339	dim_t n=pattern->numOutput;
340
341
342	colorOf=MEMALLOC(n,index_t);
343
344	if (pattern->type & PATTERN_FORMAT_SYM) {
345	Paso_setError(TYPE_ERROR,"Paso_Pattern_greedy: symmetric matrix pattern is not supported yet");
346	return;
347	}
348
349	Paso_Pattern_color(pattern,&num_colors,colorOf);
350
351	/* We do not need this loop if we set IS_IN_MIS=IS_AVAILABLE. */
352	#pragma omp parallel for private(i) schedule(static)
353	for (i=0;i<n;++i)
354	if(mis_marker[i]==IS_AVAILABLE)
355	mis_marker[i]=IS_IN_SET;
356
357	#pragma omp barrier
358	for (color=0;color<num_colors;++color) {
359	#pragma omp parallel for schedule(static) private(i,iptr,j)
360	for (i=0;i<n;++i) {
361	if (colorOf[i]==color) {
362	if (mis_marker[i]==IS_IN_SET) {
363	for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
364	j=pattern->index[iptr];
365	if (colorOf[j]<color)
366	mis_marker[j]=IS_REMOVED;
367	}
368	}
369	}
370	}
371	}
372
373
374	#pragma omp barrier
375	for (color=0;color<num_colors;++color) {
376	#pragma omp parallel for schedule(static) private(i,iptr,j)
377	for (i=0;i<n;i++) {
378	if (colorOf[i]==color) {
379	if (mis_marker[i]==IS_REMOVED) {
380	passed=TRUE;
381	for (iptr=pattern->ptr[i];iptr<pattern->ptr[i+1]; ++iptr) {
382	j=pattern->index[iptr];
383	if (colorOf[j]<color && passed) {
384	if (mis_marker[j]==IS_REMOVED) {
385	passed=TRUE;
386	}
387	else {
388	passed=FALSE;
389	/break;/
390	}
391	}
392	}
393	if (passed) mis_marker[i]=IS_IN_SET;
394	}
395
396	}
397	}
398	}
399
400	/* swap to TRUE/FALSE in mis_marker */
401	#pragma omp parallel for private(i) schedule(static)
402	for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
403
404	MEMFREE(colorOf);
405	}
406
407	#undef IS_AVAILABLE
408	#undef IS_IN_SET
409	#undef IS_REMOVED
410
411
412
413
414
415
416
417
418
419