tIPTwoBlockVec.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tIPTwoBlockVec.cc (4572B)
---
1 // Copyright (C) 2012, 2014, 2015, 2017 David Maxwell and Constantine Khroulev
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #include <cassert>
20
21 #include "IPTwoBlockVec.hh"
22 #include "pism/util/error_handling.hh"
23
24 namespace pism {
25 namespace inverse {
26
27 IPTwoBlockVec::IPTwoBlockVec(Vec a, Vec b) {
28 PetscErrorCode ierr;
29
30 MPI_Comm comm, comm_b;
31 ierr = PetscObjectGetComm((PetscObject)a, &comm); PISM_CHK(ierr, "PetscObjectGetComm");
32 ierr = PetscObjectGetComm((PetscObject)b, &comm_b); PISM_CHK(ierr, "PetscObjectGetComm");
33 assert(comm == comm_b);
34
35 // These have to be PetscInt because of PETSc calls below
36 PetscInt lo_a, hi_a;
37 ierr = VecGetOwnershipRange(a, &lo_a, &hi_a); PISM_CHK(ierr, "VecGetOwnershipRange");
38 ierr = VecGetSize(a, &m_na_global); PISM_CHK(ierr, "VecGetSize");
39 m_na_local = hi_a - lo_a;
40
41 // These have to be PetscInt because of PETSc calls below
42 PetscInt lo_b, hi_b;
43 ierr = VecGetOwnershipRange(b, &lo_b, &hi_b); PISM_CHK(ierr, "VecGetOwnershipRange");
44 ierr = VecGetSize(b, &m_nb_global); PISM_CHK(ierr, "VecGetSize");
45 m_nb_local = hi_b - lo_b;
46
47 petsc::IS is_a, is_b;
48 // a in a
49 ierr = ISCreateStride(comm, m_na_local, lo_a, 1,
50 is_a.rawptr()); PISM_CHK(ierr, "ISCreateStride");
51 // a in ab
52 ierr = ISCreateStride(comm, m_na_local, lo_a + lo_b, 1,
53 m_a_in_ab.rawptr()); PISM_CHK(ierr, "ISCreateStride");
54
55 // b in b
56 ierr = ISCreateStride(comm, m_nb_local, lo_b, 1,
57 is_b.rawptr()); PISM_CHK(ierr, "ISCreateStride");
58 // b in ab
59 ierr = ISCreateStride(comm, m_nb_local, lo_a + lo_b + m_na_local, 1,
60 m_b_in_ab.rawptr()); PISM_CHK(ierr, "ISCreateStride");
61
62 ierr = VecCreate(comm, m_ab.rawptr()); PISM_CHK(ierr, "VecCreate");
63
64 ierr = VecSetType(m_ab, "mpi"); PISM_CHK(ierr, "VecSetType");
65 ierr = VecSetSizes(m_ab, m_na_local + m_nb_local, m_na_global + m_nb_global);
66 PISM_CHK(ierr, "VecSetSizes");
67
68 ierr = VecScatterCreate(m_ab, m_a_in_ab, a, is_a,
69 m_scatter_a.rawptr()); PISM_CHK(ierr, "VecScatterCreate");
70 ierr = VecScatterCreate(m_ab, m_b_in_ab, b, is_b,
71 m_scatter_b.rawptr()); PISM_CHK(ierr, "VecScatterCreate");
72 }
73
74 IPTwoBlockVec::~IPTwoBlockVec() {
75 // empty
76 }
77
78 IS IPTwoBlockVec::blockAIndexSet() {
79 return m_a_in_ab;
80 }
81
82 IS IPTwoBlockVec::blockBIndexSet() {
83 return m_b_in_ab;
84 }
85
86 void IPTwoBlockVec::scatter(Vec a, Vec b) {
87 this->scatterToA(m_ab,a);
88 this->scatterToB(m_ab,b);
89 }
90
91 void IPTwoBlockVec::scatterToA(Vec a) {
92 this->scatterToA(m_ab,a);
93 }
94
95 void IPTwoBlockVec::scatterToB(Vec b) {
96 this->scatterToB(m_ab,b);
97 }
98
99 void IPTwoBlockVec::scatter(Vec ab, Vec a, Vec b) {
100 this->scatterToA(ab,a);
101 this->scatterToB(ab,b);
102 }
103
104 void IPTwoBlockVec::scatter_begin_end(VecScatter s, Vec a, Vec b, ScatterMode m) {
105 PetscErrorCode ierr;
106 ierr = VecScatterBegin(s, a, b, INSERT_VALUES, m);
107 PISM_CHK(ierr, "VecScatterBegin");
108
109 ierr = VecScatterEnd(s, a, b, INSERT_VALUES, m);
110 PISM_CHK(ierr, "VecScatterEnd");
111 }
112
113 void IPTwoBlockVec::scatterToA(Vec ab, Vec a) {
114 scatter_begin_end(m_scatter_a, ab, a, SCATTER_FORWARD);
115 }
116
117 void IPTwoBlockVec::scatterToB(Vec ab, Vec b) {
118 scatter_begin_end(m_scatter_b, ab, b, SCATTER_FORWARD);
119 }
120
121 void IPTwoBlockVec::gather(Vec a, Vec b) {
122 this->gatherFromA(a,m_ab);
123 this->gatherFromB(b,m_ab);
124 }
125
126 void IPTwoBlockVec::gatherFromA(Vec a) {
127 this->gatherFromA(a,m_ab);
128 }
129
130 void IPTwoBlockVec::gatherFromB(Vec b) {
131 this->gatherFromA(b,m_ab);
132 }
133
134 void IPTwoBlockVec::gather(Vec a, Vec b, Vec ab) {
135 this->gatherFromA(a,ab);
136 this->gatherFromB(b,ab);
137 }
138
139 void IPTwoBlockVec::gatherFromA(Vec a, Vec ab) {
140 scatter_begin_end(m_scatter_a, a, ab, SCATTER_REVERSE);
141 }
142
143 void IPTwoBlockVec::gatherFromB(Vec b, Vec ab) {
144 scatter_begin_end(m_scatter_b, b, ab, SCATTER_REVERSE);
145 }
146
147 } // end of namespace inverse
148 } // end of namespace pism