/[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 1917 - (show annotations)
Thu Oct 23 09:09:31 2008 UTC (11 years ago) by phornby
File MIME type: text/plain
File size: 6407 byte(s)
Remove unused vars, do a bit of standardising of indentation, and use the new Pattern_coupling.h to eliminate some implicit function declarations in Solver_AMG.c
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_MIS_NOW -2
39 #define IS_IN_MIS -3
40 #define IS_CONNECTED_TO_MIS -4
41
42
43
44 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {
45
46 index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
47 dim_t i,j;
48 double threshold=0.05;
49 index_t iptr,*index,*where_p,diagptr;
50 bool_t flag;
51 dim_t n=A->pattern->numOutput;
52 if (A->pattern->type & PATTERN_FORMAT_SYM) {
53 Paso_setError(TYPE_ERROR,"Paso_Pattern_mis: symmetric matrix pattern is not supported yet");
54 return;
55 }
56
57 /* is there any vertex available ?*/
58 while (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
59
60 #pragma omp parallel for private(i,iptr,flag) schedule(static)
61 for (i=0;i<n;++i) {
62 if (mis_marker[i]==IS_AVAILABLE) {
63 flag=IS_IN_MIS;
64 diagptr=A->pattern->ptr[i];
65 index=&(A->pattern->index[diagptr]);
66 where_p=(index_t*)bsearch(&i,
67 index,
68 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
69 sizeof(index_t),
70 Paso_comparIndex);
71 if (where_p==NULL) {
72 Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
73 } else {
74 diagptr+=(index_t)(where_p-index);
75 }
76 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
77 j=A->pattern->index[iptr]-index_offset;
78 if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {
79 flag=IS_AVAILABLE;
80 break;
81 }
82 }
83 mis_marker[i]=flag;
84 }
85 }
86
87 #pragma omp parallel for private(i,iptr) schedule(static)
88 for (i=0;i<n;i++) {
89 if (mis_marker[i]==IS_AVAILABLE) {
90 diagptr=A->pattern->ptr[i];
91 index=&(A->pattern->index[diagptr]);
92 where_p=(index_t*)bsearch(&i,
93 index,
94 A->pattern->ptr[i + 1]-A->pattern->ptr[i],
95 sizeof(index_t),
96 Paso_comparIndex);
97 if (where_p==NULL) {
98 Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
99 } else {
100 diagptr+=(index_t)(where_p-index);
101 }
102 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
103 j=A->pattern->index[iptr]-index_offset;
104 if (j!=i && mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])>=-threshold){
105 mis_marker[i]=IS_IN_MIS;
106 }
107 else {
108 mis_marker[i]=IS_CONNECTED_TO_MIS;
109 }
110 }
111 }
112 }
113 }
114 /* swap to TRUE/FALSE in mis_marker */
115 #pragma omp parallel for private(i) schedule(static)
116 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
117 }
118
119 /*
120 *
121 * Return a strength of connection mask using the classical
122 * strength of connection measure by Ruge and Stuben.
123 *
124 * Specifically, an off-diagonal entry A[i.j] is a strong
125 * connection if:
126 *
127 * -A[i,j] >= theta * max( -A[i,k] ) where k != i
128 *
129 * Otherwise, the connection is weak.
130 *
131 */
132 void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
133 {
134 index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
135 dim_t i;
136 index_t iptr;
137 double threshold,min_offdiagonal;
138 dim_t n=A->pattern->numOutput;
139 if (A->pattern->type & PATTERN_FORMAT_SYM) {
140 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
141 return;
142 }
143
144 /* is there any vertex available ?*/
145 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
146
147 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
148 for (i=0;i<n;++i) {
149 min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
150 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
151 if(A->pattern->index[iptr] != i){
152 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
153 }
154 }
155
156 threshold = theta*min_offdiagonal;
157 mis_marker[i]=IS_CONNECTED_TO_MIS;
158 #pragma omp parallel for private(iptr) schedule(static)
159 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
160 if(-1.*(A->val[iptr-index_offset]) <= threshold){
161 mis_marker[i]=IS_IN_MIS;
162 }
163 }
164 }
165 }
166 /* swap to TRUE/FALSE in mis_marker */
167 #pragma omp parallel for private(i) schedule(static)
168 for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
169 }
170 #undef IS_AVAILABLE
171 #undef IS_IN_MIS_NOW
172 #undef IS_IN_MIS
173 #undef IS_CONNECTED_TO_MIS

  ViewVC Help
Powered by ViewVC 1.1.26