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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1914 - (hide annotations)
Thu Oct 23 08:31:22 2008 UTC (10 years, 9 months ago) by phornby
File MIME type: text/plain
File size: 6488 byte(s)
Adding include files back to keep Altix happy.


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 phornby 1913 /*
32 jfenwick 1851 #include "mpi_C.h"
33     #include "Paso.h"
34     #include "Pattern.h"
35 phornby 1913 */
36 phornby 1914 #include "PasoUtil.h"
37 phornby 1913 #include "Pattern_coupling.h"
38 jfenwick 1851
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 artak 1890
48    
49 artak 1882 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker) {
50 jfenwick 1851
51     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
52 artak 1881 dim_t i,j;
53 jfenwick 1851 double threshold=0.05;
54 artak 1902 index_t iptr,*index,*where_p,diagptr;
55 jfenwick 1851 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 phornby 1907 #pragma omp parallel for private(i,iptr,flag) schedule(static)
66 jfenwick 1851 for (i=0;i<n;++i) {
67     if (mis_marker[i]==IS_AVAILABLE) {
68 artak 1902 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 jfenwick 1851 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
82 artak 1902 j=A->pattern->index[iptr]-index_offset;
83     if (j!=i && A->val[iptr]>=threshold*A->val[diagptr]) {
84     flag=IS_AVAILABLE;
85 jfenwick 1851 break;
86     }
87     }
88     mis_marker[i]=flag;
89 artak 1902 }
90     }
91    
92 phornby 1907 #pragma omp parallel for private(i,iptr) schedule(static)
93 artak 1902 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 jfenwick 1851 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
108 artak 1902 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 artak 1881 mis_marker[i]=IS_IN_MIS;
111 artak 1883 }
112     else {
113 artak 1881 mis_marker[i]=IS_CONNECTED_TO_MIS;
114 artak 1883 }
115 jfenwick 1851 }
116 artak 1902 }
117     }
118     }
119 jfenwick 1851 /* 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 artak 1890
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 phornby 1907 index_t iptr;
142 artak 1890 double threshold,min_offdiagonal;
143     bool_t flag;
144     dim_t n=A->pattern->numOutput;
145     if (A->pattern->type & PATTERN_FORMAT_SYM) {
146 artak 1902 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
147 artak 1890 return;
148     }
149    
150     /* is there any vertex available ?*/
151 artak 1902 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
152 artak 1890
153 artak 1902 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
154 artak 1890 for (i=0;i<n;++i) {
155 artak 1902 min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
156 artak 1890 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 artak 1902 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
159 artak 1890 }
160     }
161    
162     threshold = theta*min_offdiagonal;
163 artak 1902 mis_marker[i]=IS_CONNECTED_TO_MIS;
164 artak 1890 #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 artak 1902 if(-1.*(A->val[iptr-index_offset]) <= threshold){
167 artak 1890 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 jfenwick 1851 #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