/[escript]/branches/escript3047_with_pastix2995/paso/src/SystemMatrix_mergeMainAndCouple.c
ViewVC logotype

Contents of /branches/escript3047_with_pastix2995/paso/src/SystemMatrix_mergeMainAndCouple.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3049 - (show annotations)
Fri Jun 25 04:20:29 2010 UTC (8 years, 9 months ago) by lgao
File MIME type: text/plain
File size: 7709 byte(s)
Add direct solver pastix 2995. Currently works for single MPI rank only. For jobs with more than 1 MPI rank, the structure "coupleBlock" in "SystemMatrix" is not ready yet. 


1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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: SystemMatrix */
18
19 /* Merge the MainBlock and CoupleBlock in the matrix */
20
21
22 /**************************************************************/
23
24 /* Copyrights by ACcESS Australia 2003 */
25 /* Author: Lin Gao, l.gao@uq.edu.au */
26
27 /**************************************************************/
28
29 #include "Paso.h"
30 #include "SystemMatrix.h"
31
32 void Paso_SystemMatrix_mergeMainAndCouple(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val){
33 if (A->type & MATRIX_FORMAT_CSC) {
34 if (A->type & (MATRIX_FORMAT_OFFSET1 + MATRIX_FORMAT_BLK1)) {
35 Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1(A, p_ptr, p_idx, p_val);
36 } else {
37 Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_mergeMainAndCouple: CSC with index 0 or block size larger than 1 is not supported. ");
38 }
39 } else if (A->type & MATRIX_FORMAT_TRILINOS_CRS) {
40 Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_mergeMainAndCouple: TRILINOS is not supported. ");
41 } else {
42 Paso_setError(SYSTEM_ERROR,"Paso_SystemMatrix_mergeMainAndCouple: CRS is not supported. ");
43 }
44 return;
45 }
46
47 void Paso_SystemMatrix_Coupler_collectGlobalID(Paso_Coupler* coupler, index_t **global, index_t offset) {
48 Paso_MPIInfo *mpi_info = coupler->mpi_info;
49 index_t *recv_buffer=NULL;
50 index_t *send_buffer=NULL;
51 dim_t i;
52
53 if (mpi_info->size <= 1) {
54 *global = NULL;
55 return;
56 }
57
58 /* initialization */
59 send_buffer = MEMALLOC(coupler->connector->send->numSharedComponents, index_t);
60 recv_buffer = MEMALLOC(coupler->connector->recv->numSharedComponents, index_t);
61
62 /* start receiving input */
63 for (i=0; i< coupler->connector->recv->numNeighbors; ++i) {
64 #ifdef PASO_MPI
65 MPI_Irecv(&(recv_buffer[coupler->connector->recv->offsetInShared[i]]),
66 coupler->connector->recv->offsetInShared[i+1]- coupler->connector->recv->offsetInShared[i],
67 MPI_INT,
68 coupler->connector->recv->neighbor[i],
69 mpi_info->msg_tag_counter+coupler->connector->recv->neighbor[i],
70 mpi_info->comm,
71 &(coupler->mpi_requests[i]));
72 #endif
73 }
74
75 /* collect values into buffer */
76 #pragma omp parallel for private(i)
77 for (i=0; i < coupler->connector->send->numSharedComponents;++i)
78 send_buffer[i]=coupler->connector->send->shared[i] + offset;
79
80 /* send buffer out */
81 for (i=0; i< coupler->connector->send->numNeighbors; ++i) {
82 #ifdef PASO_MPI
83 MPI_Issend(&(send_buffer[coupler->connector->send->offsetInShared[i]]),
84 coupler->connector->send->offsetInShared[i+1]- coupler->connector->send->offsetInShared[i],
85 MPI_INT,
86 coupler->connector->send->neighbor[i],
87 mpi_info->msg_tag_counter+mpi_info->rank,
88 mpi_info->comm,
89 &(coupler->mpi_requests[i+ coupler->connector->recv->numNeighbors]));
90 #endif
91 }
92
93 mpi_info->msg_tag_counter+=mpi_info->size;
94
95 /* wait for receive */
96 #ifdef PASO_MPI
97 MPI_Waitall(coupler->connector->recv->numNeighbors+coupler->connector->send->numNeighbors,
98 coupler->mpi_requests,
99 coupler->mpi_stati);
100 #endif
101
102 /* finalization */
103 *global = recv_buffer;
104 MEMFREE(send_buffer);
105
106 return;
107 }
108
109 void Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val) {
110
111 index_t i, j, i_ub, j_ub, col, num_vals, main_num_vals;
112 index_t couple_num_vals, idx, rank, main_offset;
113 index_t main_num_cols=A->mainBlock->pattern->numOutput;
114 index_t couple_num_cols=A->col_coupleBlock->pattern->numOutput;
115 index_t *main_ptr=A->mainBlock->pattern->ptr;
116 index_t *main_idx=A->mainBlock->pattern->index;
117 double *main_val=A->mainBlock->val;
118 index_t *couple_ptr=A->col_coupleBlock->pattern->ptr;
119 index_t *couple_idx=A->col_coupleBlock->pattern->index;
120 double *couple_val=A->col_coupleBlock->val;
121 index_t *couple_global=NULL;
122 Paso_Coupler* coupler=A->col_coupler;
123
124 fprintf(stderr, "CHP 0\n");
125
126
127 if (A->mainBlock->col_block_size!=1 ||
128 A->mainBlock->row_block_size!=1 ||
129 A->col_coupleBlock->col_block_size!=1 ||
130 A->col_coupleBlock->row_block_size!=1) {
131 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1: requires format with block size 1.");
132 return;
133 }
134
135 if (main_num_cols != couple_num_cols) {
136 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1: number of collums do not match.");
137 return;
138 }
139
140 /* allocate arrays "ptr", "index" and "val" */
141 main_num_vals = main_ptr[main_num_cols]-1;
142 couple_num_vals = couple_ptr[couple_num_cols]-1;
143 num_vals = main_num_vals + couple_num_vals;
144 *p_ptr = MEMALLOC(main_num_cols+1, index_t);
145 *p_idx = MEMALLOC(num_vals, index_t);
146 *p_val = MEMALLOC(num_vals, double);
147
148 /* initialize before merge */
149 (*p_ptr)[0] = 1;
150 rank = A->mpi_info->rank;
151 main_offset = A->col_distribution->first_component[rank];
152 i=0;
153 j=0;
154
155 /* gather global index for entries in couple block */
156 Paso_SystemMatrix_Coupler_collectGlobalID(coupler, &couple_global, main_offset+1);
157
158 /* merge mainBlock and col_coupleBlock */
159 for (col=1; col<=main_num_cols; col++) {
160 i_ub = main_ptr[col]-1;
161 j_ub = couple_ptr[col]-1;
162 while (i < i_ub || j < j_ub) {
163 if (j < j_ub) {
164 /* switch from coupleBlock index to global row index of matrix */
165 if (A->mpi_info->size == 1) {
166 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1: requires more than 1 MPI rank when coupleBlock exists.");
167 }
168
169 idx = couple_global[couple_idx[j]-1];
170 }
171 if (j == j_ub || (i < i_ub && (main_idx[i] + main_offset) < idx)){
172 (*p_idx)[i+j] = main_idx[i] + main_offset;
173 (*p_val)[i+j] = main_val[i];
174 i++;
175 } else {
176 (*p_idx)[i+j] = idx;
177 (*p_val)[i+j] = couple_val[j];
178 j++;
179 }
180 }
181 (*p_ptr)[col] = i+j+1;
182 }
183
184 MEMFREE(couple_global);
185 return;
186 }
187
188 void Paso_SystemMatrix_copyMain_CSC_OFFSET1(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val) {
189
190 index_t i, i_ub, col, num_vals, idx;
191 index_t main_num_cols=A->mainBlock->pattern->numOutput;
192 index_t *main_ptr=A->mainBlock->pattern->ptr;
193 index_t *main_idx=A->mainBlock->pattern->index;
194 double *main_val=A->mainBlock->val;
195
196 if (A->mainBlock->col_block_size!=1 || A->mainBlock->row_block_size!=1) {
197 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1: requires format with block size 1.");
198 return;
199 }
200
201 /* allocate arrays "ptr", "index" and "val" */
202 num_vals = main_ptr[main_num_cols]-1;
203 *p_ptr = MEMALLOC(main_num_cols+1, index_t);
204 *p_idx = MEMALLOC(num_vals, index_t);
205 *p_val = MEMALLOC(num_vals, double);
206
207 /* copy from mainBlock */
208 (*p_ptr)[0] = 1;
209 i=0;
210 for (col=1; col<=main_num_cols; col++) {
211 i_ub = main_ptr[col]-1;
212 while (i < i_ub) {
213 (*p_idx)[i] = main_idx[i];
214 (*p_val)[i] = main_val[i];
215 i++;
216 }
217 (*p_ptr)[col] = i+1;
218 }
219 return;
220 }
221

  ViewVC Help
Powered by ViewVC 1.1.26