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

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

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 size: 5799 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_mergeMainAndCouple_CSC_OFFSET1(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val) {
48
49 index_t i, j, k, i_ub, j_ub, k_ub, col, num_vals, main_num_vals;
50 index_t couple_num_vals, idx, rank, tmp_idx, len;
51 index_t main_offset, couple_offset_l, couple_offset_r;
52 index_t main_num_cols=A->mainBlock->pattern->numOutput;
53 index_t couple_num_cols=A->col_coupleBlock->pattern->numOutput;
54 index_t *main_ptr=A->mainBlock->pattern->ptr;
55 index_t *main_idx=A->mainBlock->pattern->index;
56 double *main_val=A->mainBlock->val;
57 index_t *couple_ptr=A->col_coupleBlock->pattern->ptr;
58 index_t *couple_idx=A->col_coupleBlock->pattern->index;
59 double *couple_val=A->col_coupleBlock->val;
60 Paso_SharedComponents* coupler=A->col_coupler->connector->recv;
61
62 if (A->mainBlock->col_block_size!=1 ||
63 A->mainBlock->row_block_size!=1 ||
64 A->col_coupleBlock->col_block_size!=1 ||
65 A->col_coupleBlock->row_block_size!=1) {
66 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1: requires format with block size 1.");
67 return;
68 }
69
70 if (main_num_cols != couple_num_cols) {
71 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1: number of collums do not match.");
72 return;
73 }
74
75 /* allocate arrays "ptr", "index" and "val" */
76 main_num_vals = main_ptr[main_num_cols]-1;
77 couple_num_vals = couple_ptr[couple_num_cols]-1;
78 num_vals = main_num_vals + couple_num_vals;
79 *p_ptr = MEMALLOC(main_num_cols+1, index_t);
80 *p_idx = MEMALLOC(num_vals, index_t);
81 *p_val = MEMALLOC(num_vals, double);
82
83 /* initialize before merge */
84 (*p_ptr)[0] = 1;
85 rank = A->mpi_info->rank;
86 main_offset = A->col_distribution->first_component[rank];
87 k_ub=coupler->numNeighbors;
88 len = 0;
89 for (k=0; k<=k_ub; k++){
90 if (coupler->neighbor[k] < rank)
91 len += (coupler->offsetInShared[k+1] - coupler->offsetInShared[k]);
92 }
93 couple_offset_l = main_offset - len;
94 couple_offset_r = A->col_distribution->first_component[rank+1] - len;
95 i=0;
96 j=0;
97
98 /* merge mainBlock and col_coupleBlock */
99 for (col=1; col<=main_num_cols; col++) {
100 i_ub = main_ptr[col]-1;
101 j_ub = couple_ptr[col]-1;
102 while (i < i_ub || j < j_ub) {
103 if (j < j_ub) {
104 /* switch from coupleBlock index to global row index of matrix */
105 tmp_idx = coupler->shared[couple_idx[j]-1]+1-couple_num_cols;
106 for (k=1; k<=k_ub; k++)
107 if (couple_idx[j]<=coupler->offsetInShared[k]) break;
108 if (coupler->neighbor[k-1] < rank)
109 idx = couple_offset_l + tmp_idx;
110 else
111 idx = couple_offset_r + tmp_idx;
112 }
113 if (j == j_ub || (i < i_ub && (main_idx[i] + main_offset) < idx)){
114 (*p_idx)[i+j] = main_idx[i] + main_offset;
115 (*p_val)[i+j] = main_val[i];
116 i++;
117 } else {
118 (*p_idx)[i+j] = idx;
119 (*p_val)[i+j] = couple_val[j];
120 j++;
121 }
122 }
123 (*p_ptr)[col] = i+j+1;
124 }
125
126 return;
127 }
128
129 void Paso_SystemMatrix_copyMain_CSC_OFFSET1(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val) {
130
131 index_t i, i_ub, col, num_vals, idx;
132 index_t main_num_cols=A->mainBlock->pattern->numOutput;
133 index_t *main_ptr=A->mainBlock->pattern->ptr;
134 index_t *main_idx=A->mainBlock->pattern->index;
135 double *main_val=A->mainBlock->val;
136
137 if (A->mainBlock->col_block_size!=1 || A->mainBlock->row_block_size!=1) {
138 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1: requires format with block size 1.");
139 return;
140 }
141
142 /* allocate arrays "ptr", "index" and "val" */
143 num_vals = main_ptr[main_num_cols]-1;
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 /* copy from mainBlock */
149 (*p_ptr)[0] = 1;
150 i=0;
151 for (col=1; col<=main_num_cols; col++) {
152 i_ub = main_ptr[col]-1;
153 while (i < i_ub) {
154 (*p_idx)[i] = main_idx[i];
155 (*p_val)[i] = main_val[i];
156 i++;
157 }
158 (*p_ptr)[col] = i+1;
159 }
160 return;
161 }
162

  ViewVC Help
Powered by ViewVC 1.1.26