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 "mpi_C.h"
#include "Paso.h"
#include "PasoUtil.h"
#include "Pattern.h"
*/
#include "Pattern_coupling.h"


/***************************************************************/
 
#define IS_AVAILABLE -1
#define IS_IN_MIS_NOW -2
#define IS_IN_MIS -3
#define IS_CONNECTED_TO_MIS -4


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

  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 flag;
  dim_t n=A->pattern->numOutput;
  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 ?*/
     while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {

           #pragma omp parallel for private(i,iptr,flag) schedule(static) 
           for (i=0;i<n;++i) {
              if (mis_marker[i]==IS_AVAILABLE) {
                 flag=IS_IN_MIS;
                 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 && A->val[iptr]>=threshold*A->val[diagptr]) {
                        flag=IS_AVAILABLE;
                        break;
                     }
                 }
                 mis_marker[i]=flag;
                }
            }
           
              #pragma omp parallel for private(i,iptr) 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 && mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])>=-threshold){
                         mis_marker[i]=IS_IN_MIS;
                     }
                     else {
                         mis_marker[i]=IS_CONNECTED_TO_MIS;
                     }
                 }
               }
              } 
        }
     /* 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_MIS);
}

/*
 *
 * 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,j;
  index_t iptr;
  double threshold,min_offdiagonal;
  bool_t flag;
  dim_t n=A->pattern->numOutput;
  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] != i){
                min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
            }
        }

        threshold = theta*min_offdiagonal;
        mis_marker[i]=IS_CONNECTED_TO_MIS;
        #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(-1.*(A->val[iptr-index_offset]) <= threshold){
                    mis_marker[i]=IS_IN_MIS;
            }
        }
      }
    }
     /* 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_MIS);
}
#undef IS_AVAILABLE 
#undef IS_IN_MIS_NOW 
#undef IS_IN_MIS 
#undef IS_CONNECTED_TO_MIS 
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	/*
32	#include "mpi_C.h"
33	#include "Paso.h"
34	#include "PasoUtil.h"
35	#include "Pattern.h"
36	*/
37	#include "Pattern_coupling.h"
38
39
40	/***************************************************************/
41
42	#define IS_AVAILABLE -1
43	#define IS_IN_MIS_NOW -2
44	#define IS_IN_MIS -3
45	#define IS_CONNECTED_TO_MIS -4
46
47
48
49	void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {
50
51	index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
52	dim_t i,j;
53	double threshold=0.05;
54	index_t iptr,index,where_p,diagptr;
55	bool_t flag;
56	dim_t n=A->pattern->numOutput;
57	if (A->pattern->type & PATTERN_FORMAT_SYM) {
58	Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
59	return;
60	}
61
62	/* is there any vertex available ?*/
63	while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
64
65	#pragma omp parallel for private(i,iptr,flag) schedule(static)
66	for (i=0;i<n;++i) {
67	if (mis_marker[i]==IS_AVAILABLE) {
68	flag=IS_IN_MIS;
69	diagptr=A->pattern->ptr[i];
70	index=&(A->pattern->index[diagptr]);
71	where_p=(index_t*)bsearch(&i,
72	index,
73	A->pattern->ptr[i + 1]-A->pattern->ptr[i],
74	sizeof(index_t),
75	Paso_comparIndex);
76	if (where_p==NULL) {
77	Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
78	} else {
79	diagptr+=(index_t)(where_p-index);
80	}
81	for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
82	j=A->pattern->index[iptr]-index_offset;
83	if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {
84	flag=IS_AVAILABLE;
85	break;
86	}
87	}
88	mis_marker[i]=flag;
89	}
90	}
91
92	#pragma omp parallel for private(i,iptr) schedule(static)
93	for (i=0;i<n;i++) {
94	if (mis_marker[i]==IS_AVAILABLE) {
95	diagptr=A->pattern->ptr[i];
96	index=&(A->pattern->index[diagptr]);
97	where_p=(index_t*)bsearch(&i,
98	index,
99	A->pattern->ptr[i + 1]-A->pattern->ptr[i],
100	sizeof(index_t),
101	Paso_comparIndex);
102	if (where_p==NULL) {
103	Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
104	} else {
105	diagptr+=(index_t)(where_p-index);
106	}
107	for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
108	j=A->pattern->index[iptr]-index_offset;
109	if (j!=i && mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])>=-threshold){
110	mis_marker[i]=IS_IN_MIS;
111	}
112	else {
113	mis_marker[i]=IS_CONNECTED_TO_MIS;
114	}
115	}
116	}
117	}
118	}
119	/* swap to TRUE/FALSE in mis_marker */
120	#pragma omp parallel for private(i) schedule(static)
121	for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
122	}
123
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	dim_t i,j;
141	index_t iptr;
142	double threshold,min_offdiagonal;
143	bool_t flag;
144	dim_t n=A->pattern->numOutput;
145	if (A->pattern->type & PATTERN_FORMAT_SYM) {
146	Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
147	return;
148	}
149
150	/* is there any vertex available ?*/
151	if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
152
153	#pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
154	for (i=0;i<n;++i) {
155	min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
156	for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
157	if(A->pattern->index[iptr] != i){
158	min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
159	}
160	}
161
162	threshold = theta*min_offdiagonal;
163	mis_marker[i]=IS_CONNECTED_TO_MIS;
164	#pragma omp parallel for private(iptr) schedule(static)
165	for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
166	if(-1.*(A->val[iptr-index_offset]) <= threshold){
167	mis_marker[i]=IS_IN_MIS;
168	}
169	}
170	}
171	}
172	/* swap to TRUE/FALSE in mis_marker */
173	#pragma omp parallel for private(i) schedule(static)
174	for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
175	}
176	#undef IS_AVAILABLE
177	#undef IS_IN_MIS_NOW
178	#undef IS_IN_MIS
179	#undef IS_CONNECTED_TO_MIS