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"


/***************************************************************/
 
#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;*/
  index_t iptr,*index,*where_p,*diagptr;
  bool_t passed=FALSE;
  dim_t n=A->numRows;
  diagptr=MEMALLOC(n,index_t);
  /*double sum;*/
  if (A->pattern->type & PATTERN_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: 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_Solver_getAMG: 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;
            }
        }
    }
      

     /* 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);
}

/*
 *
 * Return a strength of connection mask using the classical 
 * strength of connection measure by Ruge and Stuben.
 *
 * Specifically, an off-diagonal entry A[i.j] is a strong 
 * connection if:
 *  
 *      -A[i,j] >= theta * max( -A[i,k] )   where k != i
 * 
 * Otherwise, the connection is weak.
 *  
 */  
void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta) 
{
  dim_t i;
  index_t iptr;
  double threshold,max_offdiagonal;
  dim_t 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) schedule(static) 
    for (i=0;i<n;++i) {
      if(mis_marker[i]==IS_AVAILABLE) {
        /*min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];*/
        max_offdiagonal = -1.e+15;
        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) {
            if(-(A->val[iptr])<threshold && A->pattern->index[iptr]!=i){
                    A->val[iptr]=0.;
            }
        }
       }
      }

    Paso_Pattern_coup(A,mis_marker,0.05);
}


void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta) 
{
  dim_t i,j;
  index_t iptr;
  double diag,eps_Aii,val;
  dim_t n=A->numRows;
  double* diags=MEMALLOC(n,double);
  
  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]) {
               A->val[iptr]=val;
              }
              else {
               A->val[iptr]=0.;
              }
            }
        }
       }
     }
    
    Paso_Pattern_coup(A,mis_marker,0.05);
     MEMFREE(diags);
}


#undef IS_AVAILABLE 
#undef IS_IN_SET 
#undef IS_REMOVED

void Paso_Pattern_color1(Paso_SparseMatrix* A, index_t* num_colors, index_t* colorOf) {
  index_t out=0, *mis_marker=NULL,i;
  dim_t n=A->numRows;
  mis_marker=TMPMEMALLOC(n,index_t);
  if ( !Paso_checkPtr(mis_marker) ) {
    /* get coloring */
    #pragma omp parallel for private(i) schedule(static)
    for (i = 0; i < n; ++i) {
        colorOf[i]=-1;
        mis_marker[i]=-1;
    }

    while (Paso_Util_isAny(n,colorOf,-1) && Paso_noError()) {
       Paso_Pattern_coup(A,mis_marker,0.05);

       #pragma omp parallel for private(i) schedule(static)
       for (i = 0; i < n; ++i) {
        if (mis_marker[i]) colorOf[i]=out;
        mis_marker[i]=colorOf[i];
       }
       ++out;
    }
    TMPMEMFREE(mis_marker);
    *num_colors=out;
  }
}
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
34
35	/***************************************************************/
36
37	#define IS_AVAILABLE -1
38	#define IS_IN_SET -3
39	#define IS_REMOVED -4
40
41	void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
42
43	dim_t i,j;
44	/double threshold=0.05;/
45	index_t iptr,index,where_p,*diagptr;
46	bool_t passed=FALSE;
47	dim_t n=A->numRows;
48	diagptr=MEMALLOC(n,index_t);
49	/double sum;/
50	if (A->pattern->type & PATTERN_FORMAT_SYM) {
51	Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
52	return;
53	}
54
55	#pragma omp parallel for private(i) schedule(static)
56	for (i=0;i<n;++i)
57	if(mis_marker[i]==IS_AVAILABLE)
58	mis_marker[i]=IS_IN_SET;
59
60	#pragma omp parallel for private(i,index,where_p) schedule(static)
61	for (i=0;i<n;++i) {
62	diagptr[i]=A->pattern->ptr[i];
63	index=&(A->pattern->index[diagptr[i]]);
64	where_p=(index_t*)bsearch(&i,
65	index,
66	A->pattern->ptr[i + 1]-A->pattern->ptr[i],
67	sizeof(index_t),
68	Paso_comparIndex);
69	if (where_p==NULL) {
70	Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
71	} else {
72	diagptr[i]+=(index_t)(where_p-index);
73	}
74	}
75
76
77
78	/This loop cannot be parallelized, as order matters here./
79	for (i=0;i<n;++i) {
80	if (mis_marker[i]==IS_IN_SET) {
81	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
82	j=A->pattern->index[iptr];
83	if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr[i]])) {
84	mis_marker[j]=IS_REMOVED;
85	}
86	}
87	}
88	}
89
90
91
92	/This loop cannot be parallelized, as order matters here./
93	for (i=0;i<n;i++) {
94	if (mis_marker[i]==IS_REMOVED) {
95	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
96	j=A->pattern->index[iptr];
97	if (mis_marker[j]==IS_IN_SET) {
98	if ((A->val[iptr]/A->val[diagptr[i]])>=-threshold) {
99	passed=TRUE;
100	}
101	else {
102	passed=FALSE;
103	break;
104	}
105	}
106	}
107	if (passed) {
108	mis_marker[i]=IS_IN_SET;
109	passed=FALSE;
110	}
111	}
112	}
113
114
115	/* swap to TRUE/FALSE in mis_marker */
116	#pragma omp parallel for private(i) schedule(static)
117	for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]!=IS_IN_SET);
118
119	MEMFREE(diagptr);
120	}
121
122	/*
123	*
124	* Return a strength of connection mask using the classical
125	* strength of connection measure by Ruge and Stuben.
126	*
127	* Specifically, an off-diagonal entry A[i.j] is a strong
128	* connection if:
129	*
130	* -A[i,j] >= theta * max( -A[i,k] ) where k != i
131	*
132	* Otherwise, the connection is weak.
133	*
134	*/
135	void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
136	{
137	dim_t i;
138	index_t iptr;
139	double threshold,max_offdiagonal;
140	dim_t n=A->numRows;
141	if (A->pattern->type & PATTERN_FORMAT_SYM) {
142	Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
143	return;
144	}
145
146	#pragma omp parallel for private(i,iptr,max_offdiagonal,threshold) schedule(static)
147	for (i=0;i<n;++i) {
148	if(mis_marker[i]==IS_AVAILABLE) {
149	/min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];/
150	max_offdiagonal = -1.e+15;
151	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
152	if(A->pattern->index[iptr] != i){
153	max_offdiagonal = MAX(max_offdiagonal,-(A->val[iptr]));
154	}
155	}
156
157	threshold = theta*max_offdiagonal;
158	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
159	if(-(A->val[iptr])<threshold && A->pattern->index[iptr]!=i){
160	A->val[iptr]=0.;
161	}
162	}
163	}
164	}
165
166	Paso_Pattern_coup(A,mis_marker,0.05);
167	}
168
169
170
171
172	void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
173	{
174	dim_t i,j;
175	index_t iptr;
176	double diag,eps_Aii,val;
177	dim_t n=A->numRows;
178	double* diags=MEMALLOC(n,double);
179
180	if (A->pattern->type & PATTERN_FORMAT_SYM) {
181	Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
182	return;
183	}
184
185
186	#pragma omp parallel for private(i,iptr,diag) schedule(static)
187	for (i=0;i<n;++i) {
188	diag = 0;
189	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
190	if(A->pattern->index[iptr] != i){
191	diag+=A->val[iptr];
192	}
193	}
194	diags[i]=ABS(diag);
195	}
196
197
198	#pragma omp parallel for private(i,iptr,j,val,eps_Aii) schedule(static)
199	for (i=0;i<n;++i) {
200	if (mis_marker[i]==IS_AVAILABLE) {
201	eps_Aii = thetathetadiags[i];
202	val=0.;
203	for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
204	j=A->pattern->index[iptr];
205	val=A->val[iptr];
206	if(j!= i) {
207	if(valval>=eps_Aii diags[j]) {
208	A->val[iptr]=val;
209	}
210	else {
211	A->val[iptr]=0.;
212	}
213	}
214	}
215	}
216	}
217
218	Paso_Pattern_coup(A,mis_marker,0.05);
219	MEMFREE(diags);
220	}
221
222
223	#undef IS_AVAILABLE
224	#undef IS_IN_SET
225	#undef IS_REMOVED
226
227	void Paso_Pattern_color1(Paso_SparseMatrix* A, index_t* num_colors, index_t* colorOf) {
228	index_t out=0, *mis_marker=NULL,i;
229	dim_t n=A->numRows;
230	mis_marker=TMPMEMALLOC(n,index_t);
231	if ( !Paso_checkPtr(mis_marker) ) {
232	/* get coloring */
233	#pragma omp parallel for private(i) schedule(static)
234	for (i = 0; i < n; ++i) {
235	colorOf[i]=-1;
236	mis_marker[i]=-1;
237	}
238
239	while (Paso_Util_isAny(n,colorOf,-1) && Paso_noError()) {
240	Paso_Pattern_coup(A,mis_marker,0.05);
241
242	#pragma omp parallel for private(i) schedule(static)
243	for (i = 0; i < n; ++i) {
244	if (mis_marker[i]) colorOf[i]=out;
245	mis_marker[i]=colorOf[i];
246	}
247	++out;
248	}
249	TMPMEMFREE(mis_marker);
250	*num_colors=out;
251	}
252	}