/[escript]/branches/doubleplusgood/paso/src/SystemMatrix_debug.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/SystemMatrix_debug.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4291 - (show annotations)
Thu Mar 7 08:57:58 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 6483 byte(s)


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au
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 Paso: SystemMatrix: debugging tools
16
17 ************************************************************************************
18
19 Author: Lutz Gross, l.gross@uq.edu.au
20
21 ************************************************************************************/
22
23 #include "SystemMatrix.h"
24 #include "esysUtils/error.h"
25
26 /************************************************************************************/
27
28 /* fills the matrix with values i+f1*j where i and j are the global row
29 * and column indices of the matrix entry */
30 void Paso_SystemMatrix_fillWithGlobalCoordinates(Paso_SystemMatrix *A, const double f1)
31 {
32 dim_t ib, iPtr, q, i,p;
33 const dim_t n=Paso_SystemMatrix_getNumRows(A);
34 const dim_t m=Paso_SystemMatrix_getNumCols(A);
35 const index_t me = A->mpi_info->rank;
36 const dim_t block_size=A->block_size;
37 const index_t row_offset = A->row_distribution->first_component[me];
38 const index_t col_offset = A->col_distribution->first_component[me];
39 Paso_Coupler* col_coupler=NULL, *row_coupler=NULL;
40 double *cols=NULL, *rows=NULL;
41
42 cols=new double[m];
43 rows=new double[n];
44 col_coupler= Paso_Coupler_alloc(A->col_coupler->connector, 1);
45 row_coupler = Paso_Coupler_alloc(A->col_coupler->connector, 1);
46
47 #pragma omp parallel for private(i)
48 for (i=0; i<n; ++i) rows[i]=row_offset+i;
49 Paso_Coupler_startCollect(col_coupler, rows);
50
51 #pragma omp parallel for private(i)
52 for (i=0; i<m; ++i) cols[i]=col_offset+i;
53 Paso_Coupler_startCollect(row_coupler, cols);
54
55
56
57 /* main block : */
58 for (q=0; q< n; ++q){
59 for (iPtr =A->mainBlock->pattern->ptr[q]; iPtr<A->mainBlock->pattern->ptr[q+1]; ++iPtr) {
60 p =A->mainBlock->pattern->index[iPtr];
61 for (ib=0; ib<block_size; ib++) A->mainBlock->val[iPtr*block_size+ib]=f1*rows[q]+cols[p];
62 }
63 }
64
65 Paso_Coupler_finishCollect(col_coupler);
66 if ( A->col_coupleBlock != NULL ){
67 for (q=0; q< A->col_coupleBlock->pattern->numOutput; ++q){
68 for (iPtr =A->col_coupleBlock->pattern->ptr[q]; iPtr<A->col_coupleBlock->pattern->ptr[q+1]; ++iPtr) {
69 p =A->col_coupleBlock->pattern->index[iPtr];
70 for (ib=0; ib<block_size; ib++) A->col_coupleBlock->val[iPtr*block_size+ib]=f1*rows[q]+col_coupler->recv_buffer[p];
71 }
72 }
73 }
74
75 Paso_Coupler_finishCollect(row_coupler);
76 if ( A->row_coupleBlock != NULL ){
77 for (p=0; p< A->row_coupleBlock->pattern->numOutput; ++p){
78 for (iPtr =A->row_coupleBlock->pattern->ptr[p]; iPtr<A->row_coupleBlock->pattern->ptr[p+1]; ++iPtr) {
79 q =A->row_coupleBlock->pattern->index[iPtr];
80 for (ib=0; ib<block_size; ib++) A->row_coupleBlock->val[iPtr*block_size+ib]=row_coupler->recv_buffer[p]*f1+cols[q];
81 }
82 }
83 }
84
85 delete[] cols;
86 delete[] rows;
87 Paso_Coupler_free(row_coupler);
88 Paso_Coupler_free(col_coupler);
89 return;
90 }
91
92 void Paso_SystemMatrix_print(Paso_SystemMatrix *A)
93 {
94 dim_t iPtr, q, p, ib;
95 const dim_t n=Paso_SystemMatrix_getNumRows(A);
96 const dim_t block_size=A->block_size;
97 index_t rank=A->mpi_info->rank;
98 char *str1, *str2;
99 str1 = new char[n*n*block_size*30+100];
100 str2 = new char[30];
101
102 sprintf(str1, "rank %d Main Block:\n-----------\n", rank);
103 for (q=0; q< n; ++q){
104 sprintf(str2, "Row %d: ",q);
105 strcat(str1, str2);
106 for (iPtr =A->mainBlock->pattern->ptr[q]; iPtr<A->mainBlock->pattern->ptr[q+1]; ++iPtr) {
107 sprintf(str2, "(%d ",A->mainBlock->pattern->index[iPtr]);
108 strcat(str1, str2);
109 for (ib=0; ib<block_size; ib++){
110 sprintf(str2, "%f ", A->mainBlock->val[iPtr*block_size+ib]);
111 strcat(str1, str2);
112 }
113 sprintf(str2, "),");
114 strcat(str1, str2);
115 }
116 sprintf(str1, "%s\n", str1);
117 }
118 fprintf(stderr, "%s", str1);
119 if ( ( A->col_coupleBlock != NULL) && (A->mpi_info->size>1) ) {
120 sprintf(str1, "rank %d Column Couple Block:\n------------------\n", rank);
121 for (q=0; q< A->col_coupleBlock->pattern->numOutput; ++q){
122 sprintf(str2, "Row %d: ", q);
123 strcat(str1, str2);
124 for (iPtr =A->col_coupleBlock->pattern->ptr[q]; iPtr<A->col_coupleBlock->pattern->ptr[q+1]; ++iPtr) {
125 if (A->global_id)
126 sprintf(str2, "(%d %f),",A->global_id[A->col_coupleBlock->pattern->index[iPtr]], A->col_coupleBlock->val[iPtr*block_size]);
127 else
128 sprintf(str2, "(%d %f),",A->col_coupleBlock->pattern->index[iPtr], A->col_coupleBlock->val[iPtr*block_size]);
129 strcat(str1, str2);
130 }
131 sprintf(str1, "%s\n", str1);
132 }
133 fprintf(stderr, "%s", str1);
134 }
135 if ( ( A->row_coupleBlock != NULL ) && (A->mpi_info->size>1) ) {
136 sprintf(str1, "rank %d Row Couple Block:\n--------------------\n", rank);
137 for (p=0; p< A->row_coupleBlock->pattern->numOutput; ++p){
138 sprintf(str2, "Row %d:", p);
139 strcat(str1, str2);
140 for (iPtr =A->row_coupleBlock->pattern->ptr[p]; iPtr<A->row_coupleBlock->pattern->ptr[p+1]; ++iPtr) {
141 sprintf(str2, "(%d %f),",A->row_coupleBlock->pattern->index[iPtr], A->row_coupleBlock->val[iPtr*block_size]);
142 strcat(str1, str2);
143 }
144 sprintf(str1, "%s\n", str1);
145 }
146 fprintf(stderr, "%s", str1);
147 }
148 if ( ( A->remote_coupleBlock != NULL ) && (A->mpi_info->size>1) ) {
149 sprintf(str1, "rank %d Remote Couple Block:\n--------------------\n", rank);
150 for (p=0; p< A->remote_coupleBlock->pattern->numOutput; ++p){
151 sprintf(str2, "Row %d:", p);
152 strcat(str1, str2);
153 for (iPtr =A->remote_coupleBlock->pattern->ptr[p]; iPtr<A->remote_coupleBlock->pattern->ptr[p+1]; ++iPtr) {
154 sprintf(str2, "(%d %f),",A->remote_coupleBlock->pattern->index[iPtr], A->remote_coupleBlock->val[iPtr*block_size]);
155 strcat(str1, str2);
156 }
157 sprintf(str1, "%s\n", str1);
158 }
159 fprintf(stderr, "%s", str1);
160 }
161 delete[] str1;
162 delete[] str2;
163 return;
164 }
165

Properties

Name Value
svn:executable *

  ViewVC Help
Powered by ViewVC 1.1.26