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 "Solver.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 naib,iptr;
  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(naib,i,iptr,flag) schedule(static) 
           for (i=0;i<n;++i) {
              if (mis_marker[i]==IS_AVAILABLE) {
                 flag=IS_AVAILABLE;
                 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
                     naib=A->pattern->index[iptr]-index_offset;
                     if (naib!=i && A->val[iptr]>=threshold*A->val[i]) {
                        flag=IS_IN_MIS;
                        break;
                     }
                 }
                 mis_marker[i]=flag;
              }
           }

           #pragma omp parallel for private(naib,i,iptr) schedule(static)
           for (i=0;i<n;i++) {
              if (mis_marker[i]==IS_AVAILABLE) {
                 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
                     naib=A->pattern->index[iptr]-index_offset;
                     if (naib!=i && mis_marker[naib]==IS_IN_MIS && A->val[iptr]/A->val[i]>=-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 naib,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_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,min_offdiagonal,threshold) schedule(static) 
      for (i=0;i<n;++i) {
        min_offdiagonal = A->val[A->pattern->ptr[i]];
        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]);
            }
        }

        threshold = theta*min_offdiagonal;
        #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) {
            mis_marker[i]=IS_CONNECTED_TO_MIS;
            if(A->val[iptr] <= threshold){
                if(A->pattern->index[iptr] != i){
                    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	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			#include "mpi_C.h"
32			#include "Paso.h"
33			#include "PasoUtil.h"
34			#include "Pattern.h"
35			#include "Solver.h"
36
37
38			/***************************************************************/
39
40			#define IS_AVAILABLE -1
41			#define IS_IN_MIS_NOW -2
42			#define IS_IN_MIS -3
43			#define IS_CONNECTED_TO_MIS -4
44
45	artak	1890
46
47	artak	1882	void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {
48	jfenwick	1851
49			index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
50	artak	1881	dim_t i,j;
51	jfenwick	1851	double threshold=0.05;
52			index_t naib,iptr;
53			bool_t flag;
54			dim_t n=A->pattern->numOutput;
55			if (A->pattern->type & PATTERN_FORMAT_SYM) {
56			Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
57			return;
58			}
59
60			/* is there any vertex available ?*/
61			while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
62
63			#pragma omp parallel for private(naib,i,iptr,flag) schedule(static)
64			for (i=0;i<n;++i) {
65			if (mis_marker[i]==IS_AVAILABLE) {
66	artak	1881	flag=IS_AVAILABLE;
67	jfenwick	1851	for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
68			naib=A->pattern->index[iptr]-index_offset;
69	artak	1883	if (naib!=i && A->val[iptr]>=threshold*A->val[i]) {
70	artak	1881	flag=IS_IN_MIS;
71	jfenwick	1851	break;
72			}
73			}
74			mis_marker[i]=flag;
75			}
76			}
77
78			#pragma omp parallel for private(naib,i,iptr) schedule(static)
79			for (i=0;i<n;i++) {
80			if (mis_marker[i]==IS_AVAILABLE) {
81			for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
82			naib=A->pattern->index[iptr]-index_offset;
83	artak	1883	if (naib!=i && mis_marker[naib]==IS_IN_MIS && A->val[iptr]/A->val[i]>=-threshold){
84	artak	1881	mis_marker[i]=IS_IN_MIS;
85	artak	1883	}
86			else {
87	artak	1881	mis_marker[i]=IS_CONNECTED_TO_MIS;
88	artak	1883	}
89	jfenwick	1851	}
90			}
91			}
92			}
93			/* swap to TRUE/FALSE in mis_marker */
94			#pragma omp parallel for private(i) schedule(static)
95			for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
96			}
97	artak	1890
98			/*
99			*
100			* Return a strength of connection mask using the classical
101			* strength of connection measure by Ruge and Stuben.
102			*
103			* Specifically, an off-diagonal entry A[i.j] is a strong
104			* connection if:
105			*
106			* -A[i,j] >= theta * max( -A[i,k] ) where k != i
107			*
108			* Otherwise, the connection is weak.
109			*
110			*/
111			void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
112			{
113			index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
114			dim_t i,j;
115			index_t naib,iptr;
116			double threshold,min_offdiagonal;
117			bool_t flag;
118			dim_t n=A->pattern->numOutput;
119			if (A->pattern->type & PATTERN_FORMAT_SYM) {
120			Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
121			return;
122			}
123
124			/* is there any vertex available ?*/
125			while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
126
127			#pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
128			for (i=0;i<n;++i) {
129			min_offdiagonal = A->val[A->pattern->ptr[i]];
130			for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
131			if(A->pattern->index[iptr] != i){
132			min_offdiagonal = MIN(min_offdiagonal,A->val[iptr]);
133			}
134			}
135
136			threshold = theta*min_offdiagonal;
137			#pragma omp parallel for private(iptr) schedule(static)
138			for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
139			mis_marker[i]=IS_CONNECTED_TO_MIS;
140			if(A->val[iptr] <= threshold){
141			if(A->pattern->index[iptr] != i){
142			mis_marker[i]=IS_IN_MIS;
143			}
144			}
145			}
146			}
147			}
148			/* swap to TRUE/FALSE in mis_marker */
149			#pragma omp parallel for private(i) schedule(static)
150			for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
151			}
152	jfenwick	1851	#undef IS_AVAILABLE
153			#undef IS_IN_MIS_NOW
154			#undef IS_IN_MIS
155			#undef IS_CONNECTED_TO_MIS