/[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 1913 - (show annotations)
Thu Oct 23 08:27:33 2008 UTC (10 years, 9 months ago) by phornby
File MIME type: text/plain
File size: 6488 byte(s)
Create Pattern_coupling.h as there were too many modules that implicitly defined these
functions. This checkin is to test the minimum number of files this change allows me to
include to get Pattern_coupling.c to compile correctly on the Altix.


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

  ViewVC Help
Powered by ViewVC 1.1.26