/[escript]/trunk/paso/src/Pattern_coupling.c
ViewVC logotype

Contents of /trunk/paso/src/Pattern_coupling.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1890 - (show annotations)
Fri Oct 17 00:14:22 2008 UTC (11 years ago) by artak
File MIME type: text/plain
File size: 5214 byte(s)
Current version of AMG. check levels for stoping applies presmoothing, but not ready yet
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 "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
46
47 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {
48
49 index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
50 dim_t i,j;
51 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 flag=IS_AVAILABLE;
67 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 if (naib!=i && A->val[iptr]>=threshold*A->val[i]) {
70 flag=IS_IN_MIS;
71 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 if (naib!=i && mis_marker[naib]==IS_IN_MIS && A->val[iptr]/A->val[i]>=-threshold){
84 mis_marker[i]=IS_IN_MIS;
85 }
86 else {
87 mis_marker[i]=IS_CONNECTED_TO_MIS;
88 }
89 }
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
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 #undef IS_AVAILABLE
153 #undef IS_IN_MIS_NOW
154 #undef IS_IN_MIS
155 #undef IS_CONNECTED_TO_MIS

  ViewVC Help
Powered by ViewVC 1.1.26