/[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 1931 - (hide annotations)
Mon Oct 27 01:31:28 2008 UTC (11 years, 2 months ago) by artak
File MIME type: text/plain
File size: 6634 byte(s)
Direct solver is added to solve_AMG. It is now working but I guess need optimization.
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 1914 #include "PasoUtil.h"
32 phornby 1913 #include "Pattern_coupling.h"
33 jfenwick 1851
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 artak 1890
43    
44 artak 1931 void Paso_Pattern_coup(Paso_SparseMatrix* A, index_t* mis_marker, double threshold) {
45 jfenwick 1851
46     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
47 artak 1881 dim_t i,j;
48 artak 1931 /*double threshold=0.05;*/
49 artak 1902 index_t iptr,*index,*where_p,diagptr;
50 artak 1931 bool_t flag,fail;
51 jfenwick 1851 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 artak 1931 #pragma omp parallel for private(i,iptr,index,where_p,diagptr,j) schedule(static)
61     for (i=0;i<n;++i) {
62 jfenwick 1851 if (mis_marker[i]==IS_AVAILABLE) {
63 artak 1902 diagptr=A->pattern->ptr[i];
64     index=&(A->pattern->index[diagptr]);
65     where_p=(index_t*)bsearch(&i,
66     index,
67     A->pattern->ptr[i + 1]-A->pattern->ptr[i],
68     sizeof(index_t),
69     Paso_comparIndex);
70     if (where_p==NULL) {
71     Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
72     } else {
73     diagptr+=(index_t)(where_p-index);
74     }
75 jfenwick 1851 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
76 artak 1902 j=A->pattern->index[iptr]-index_offset;
77 artak 1931 if (j!=i && ABS(A->val[iptr])>=threshold*ABS(A->val[diagptr])) {
78     mis_marker[j]=IS_CONNECTED_TO_MIS;
79 jfenwick 1851 }
80     }
81 artak 1902 }
82     }
83 artak 1931
84     #pragma omp parallel for private(i) schedule(static)
85     for (i=0;i<n;++i)
86     if(mis_marker[i]==IS_AVAILABLE)
87     mis_marker[i]=IS_IN_MIS;
88 artak 1902
89 artak 1931 #pragma omp parallel for private(i,iptr,fail,index,where_p,diagptr) schedule(static)
90 artak 1902 for (i=0;i<n;i++) {
91 artak 1931 if (mis_marker[i]==IS_CONNECTED_TO_MIS) {
92     fail=FALSE;
93 artak 1902 diagptr=A->pattern->ptr[i];
94     index=&(A->pattern->index[diagptr]);
95     where_p=(index_t*)bsearch(&i,
96     index,
97     A->pattern->ptr[i + 1]-A->pattern->ptr[i],
98     sizeof(index_t),
99     Paso_comparIndex);
100     if (where_p==NULL) {
101     Paso_setError(VALUE_ERROR, "Paso_Solver_getAMG: main diagonal element missing.");
102     } else {
103     diagptr+=(index_t)(where_p-index);
104     }
105 jfenwick 1851 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
106 artak 1902 j=A->pattern->index[iptr]-index_offset;
107 artak 1931 if (mis_marker[j]==IS_IN_MIS && (A->val[iptr]/A->val[diagptr])<-threshold){
108     fail=TRUE;
109     break;
110 artak 1883 }
111 jfenwick 1851 }
112 artak 1931 if(!fail)
113     mis_marker[i]=IS_IN_MIS;
114 artak 1902 }
115     }
116     }
117 jfenwick 1851 /* swap to TRUE/FALSE in mis_marker */
118     #pragma omp parallel for private(i) schedule(static)
119     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
120     }
121 artak 1890
122     /*
123     *
124     * Return a strength of connection mask using the classical
125     * strength of connection measure by Ruge and Stuben.
126     *
127     * Specifically, an off-diagonal entry A[i.j] is a strong
128     * connection if:
129     *
130     * -A[i,j] >= theta * max( -A[i,k] ) where k != i
131     *
132     * Otherwise, the connection is weak.
133     *
134     */
135     void Paso_Pattern_RS(Paso_SparseMatrix* A, index_t* mis_marker, double theta)
136     {
137     index_t index_offset=(A->pattern->type & PATTERN_FORMAT_OFFSET1 ? 1:0);
138 phornby 1917 dim_t i;
139 phornby 1907 index_t iptr;
140 artak 1890 double threshold,min_offdiagonal;
141     dim_t n=A->pattern->numOutput;
142     if (A->pattern->type & PATTERN_FORMAT_SYM) {
143 artak 1902 Paso_setError(TYPE_ERROR,"Paso_Pattern_RS: symmetric matrix pattern is not supported yet");
144 artak 1890 return;
145     }
146    
147     /* is there any vertex available ?*/
148 artak 1902 if (Paso_Util_isAny(n,mis_marker,IS_AVAILABLE)) {
149 artak 1890
150 artak 1902 #pragma omp parallel for private(i,iptr,min_offdiagonal,threshold) schedule(static)
151 artak 1890 for (i=0;i<n;++i) {
152 artak 1902 min_offdiagonal = A->val[A->pattern->ptr[i]-index_offset];
153 artak 1890 for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
154     if(A->pattern->index[iptr] != i){
155 artak 1902 min_offdiagonal = MIN(min_offdiagonal,A->val[iptr-index_offset]);
156 artak 1890 }
157     }
158    
159     threshold = theta*min_offdiagonal;
160 artak 1931 mis_marker[i]=IS_IN_MIS;
161 artak 1890 #pragma omp parallel for private(iptr) schedule(static)
162     for (iptr=A->pattern->ptr[i]-index_offset;iptr<A->pattern->ptr[i+1]-index_offset; ++iptr) {
163 artak 1931 if(A->val[iptr-index_offset] < threshold){
164     mis_marker[i]=IS_CONNECTED_TO_MIS;
165     break;
166 artak 1890 }
167     }
168     }
169     }
170     /* swap to TRUE/FALSE in mis_marker */
171     #pragma omp parallel for private(i) schedule(static)
172     for (i=0;i<n;i++) mis_marker[i]=(mis_marker[i]==IS_IN_MIS);
173     }
174 jfenwick 1851 #undef IS_AVAILABLE
175     #undef IS_IN_MIS_NOW
176     #undef IS_IN_MIS
177     #undef IS_CONNECTED_TO_MIS

  ViewVC Help
Powered by ViewVC 1.1.26