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) {

  index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
  dim_t i,j;
  /*double threshold=0.05;*/
  index_t iptr,*index,*where_p,diagptr;
  bool_t fail;
  dim_t n=A->numRows;
  /*double sum;*/
  if (A->pattern->type & PATTERN_FORMAT_SYM) {
    Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
    return;
  }
   
     /* is there any vertex available ?*/
     if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {

            #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static) 
            for (i=0;i<n;++i) {
              if (mis_marker[i]==IS_AVAILABLE) {
                 diagptr=A->pattern->ptr[i];
                 index=&(A->pattern->index[diagptr]);
                 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+=(index_t)(where_p-index);
                }
                 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
                     j=A->pattern->index[iptr]-index_offset;
                     if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
                        mis_marker[j]=IS_REMOVED;
                       }
                 }
                }
            }
            
            #pragma omp parallel for private(i,iptr) 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,iptr,index,where_p,diagptr) schedule(static)
              for (i=0;i<n;i++) {
               if (mis_marker[i]==IS_REMOVED) {
                 fail=FALSE;
                 diagptr=A->pattern->ptr[i];
                 index=&(A->pattern->index[diagptr]);
                 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+=(index_t)(where_p-index);
                }
                 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
                     j=A->pattern->index[iptr]-index_offset;
                     if (mis_marker[j]==IS_IN_SET && (A->val[iptr]/A->val[diagptr])<-threshold){
                         fail=TRUE;
                         #ifndef _OPENMP  
                         break;
                         #endif
                     }
                 }
                 if(!fail) {
                    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);
}

/*
 *
 * 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) 
{
  index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
  dim_t i;
  index_t iptr;
  double threshold,min_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;
  }

/* is there any vertex available ?*/
     if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {

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

        threshold = theta*min_offdiagonal;
        mis_marker[i]=IS_IN_SET;
        #pragma omp parallel for private(iptr) schedule(static) 
        for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
            if(A->val[iptr-index_offset] < threshold){
                    mis_marker[i]=IS_REMOVED;
                    #ifndef _OPENMP    
                     break;
                    #endif
            }
        }
      }
    }
     /* 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);
}


void Paso_Pattern_Aggregiation(Paso_SparseMatrix* A, index_t* mis_marker, double theta) 
{
  index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
  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;
  }
    
   if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
    #pragma omp parallel for private(i,iptr,diag) schedule(static) 
      for (i=0;i<n;++i) {
        diag = 0;
        for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
            if(A->pattern->index[iptr-index_offset] != i){
                diag+=A->val[iptr-index_offset];
            }
        }
        diags[i]=ABS(diag);
      }
    
    #pragma omp parallel for private(i,iptr,diag,j,val,eps_Aii) schedule(static) 
      for (i=0;i<n;++i) {
        eps_Aii = theta*theta*diags[i];
        mis_marker[i]=IS_REMOVED;

        for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
            j=A->pattern->index[iptr-index_offset];
            val=A->val[iptr-index_offset];
            if(j!= i && val*val>=eps_Aii * diags[j]){
                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(diags);
}


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