PetscSolverFetiOperators.cc 15.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.

#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiStructs.h"
#include "parallel/PetscSolverFetiTimings.h"

namespace AMDiS {

  int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
  {
    // S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi

    void *ctx;
    MatShellGetContext(mat, &ctx);
    SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);

    MatMult(data->subSolver->getMatInteriorCoarse(), x, data->tmp_vec_b);
    data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
    MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b, 
	    data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarse(), x, y);
    VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);

    return 0;
  }


37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
  int petscMultMatSchurPrimalAugmented(Mat mat, Vec x, Vec y)
  {
    void *ctx;
    MatShellGetContext(mat, &ctx);
    SchurPrimalAugmentedData* data = 
      static_cast<SchurPrimalAugmentedData*>(ctx);

    Vec x_primal, x_mu, y_primal, y_mu;    
    VecNestGetSubVec(x, 0, &x_primal);
    VecNestGetSubVec(x, 1, &x_mu);
    VecNestGetSubVec(y, 0, &y_primal);
    VecNestGetSubVec(y, 1, &y_mu);

    // inv(K_BB) K_BPi x_Pi
    MatMult(data->subSolver->getMatInteriorCoarse(), x_primal, data->tmp_vec_b0);
    data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);

    // inv(K_BB) J^T Q x_mu
    MatMult(*(data->mat_augmented_lagrange), x_mu, data->tmp_vec_lagrange);
    MatMultTranspose(*(data->mat_lagrange), data->tmp_vec_lagrange, data->tmp_vec_b1);
    data->subSolver->solveGlobal(data->tmp_vec_b1, data->tmp_vec_b1);

    // y_Pi = (K_PiPi - K_PiB inv(K_BB) K_BPi) x_pi
    MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b0,
	    data->tmp_vec_primal);
    MatMult(data->subSolver->getMatCoarse(), x_primal, y_primal);
    VecAXPY(y_primal, -1.0, data->tmp_vec_primal);

    // y_Pi += (-K_PiB inv(K_BB) J^T Q) x_mu
    MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b1,
	    data->tmp_vec_primal);
    VecAXPY(y_primal, -1.0, data->tmp_vec_primal);

    // y_mu = (-Q^T J inv(K_BB) K_BPi) x_Pi + (-Q^T J inv(K_BB) J^T Q) x_mu
    //      = -Q^T J (inv(K_BB) K_BPi x_Pi + inv(K_BB) J^T Q x_mu)
    VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
74
    MatMultTranspose(*(data->mat_augmented_lagrange), data->tmp_vec_lagrange, y_mu);
75
76
77
78
79
80
    VecScale(y_mu, -1.0);

    return 0;
  }


81
82
83
84
85
  // y = mat * x
  int petscMultMatFeti(Mat mat, Vec x, Vec y)
  {
    FUNCNAME("petscMultMatFeti()");

86
87
    //    F = J inv(K_BB) trans(J) + J inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(J)
    // => F = J [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(J)
88
89
90
91
92
93
94

    double wtime = MPI::Wtime();

    void *ctx;
    MatShellGetContext(mat, &ctx);
    FetiData* data = static_cast<FetiData*>(ctx);

Thomas Witkowski's avatar
Thomas Witkowski committed
95
    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
96
97

    double wtime01 = MPI::Wtime();
Thomas Witkowski's avatar
Thomas Witkowski committed
98
    data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
99
100
101

    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
102
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
103
104

    MatMult(data->subSolver->getMatCoarseInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
105
	    data->tmp_vec_b0, data->tmp_vec_primal0);
106
107

    wtime01 = MPI::Wtime();
Thomas Witkowski's avatar
Thomas Witkowski committed
108
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal0, data->tmp_vec_primal0);
109
110
111
    FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);

    MatMult(data->subSolver->getMatInteriorCoarse(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
112
	    data->tmp_vec_primal0, data->tmp_vec_b0);
113
114

    wtime01 = MPI::Wtime();
Thomas Witkowski's avatar
Thomas Witkowski committed
115
    data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
116
117
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
118
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y);
119
120
121
122
123
124
125
126
127

    VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);

    FetiTimings::fetiSolve += (MPI::Wtime() - wtime);

    return 0;
  }


128
129
130
131
132
133
134
135
136
137
  // y = mat * x
  int petscMultMatFetiAugmented(Mat mat, Vec x, Vec y)
  {
    FUNCNAME("petscMultMatFetiAugmented()");

    void *ctx;
    MatShellGetContext(mat, &ctx);
    FetiData* data = static_cast<FetiData*>(ctx);

    Vec vec_mu0, vec_mu1;
Thomas Witkowski's avatar
Thomas Witkowski committed
138
    MatGetVecs(*(data->mat_augmented_lagrange), &vec_mu0, PETSC_NULL);
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
    VecDuplicate(vec_mu0, &vec_mu1);

    MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
    data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);

    MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b0, data->tmp_vec_primal0);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
    MatMultTranspose(*(data->mat_augmented_lagrange), data->tmp_vec_lagrange, vec_mu0);

    Vec vec_array0[2] = {data->tmp_vec_primal0, vec_mu0};
    Vec vec_array1[2] = {data->tmp_vec_primal1, vec_mu1};
    Vec vec_nest0, vec_nest1;
    VecCreateNest(PETSC_COMM_WORLD, 2, PETSC_NULL, vec_array0, &vec_nest0);
    VecCreateNest(PETSC_COMM_WORLD, 2, PETSC_NULL, vec_array1, &vec_nest1);

    KSPSolve(*(data->ksp_schur_primal), vec_nest0, vec_nest1);

    // Step 1
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y);

    // Step 2
    MatMult(data->subSolver->getMatInteriorCoarse(), data->tmp_vec_primal1, data->tmp_vec_b0);
    data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
    VecAXPY(y, 1.0, data->tmp_vec_lagrange);

    // Step 3
    MatMult(*(data->mat_augmented_lagrange), vec_mu1, data->tmp_vec_lagrange);
    MatMultTranspose(*(data->mat_lagrange), data->tmp_vec_lagrange, data->tmp_vec_b0);
    data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
    VecAXPY(y, 1.0, data->tmp_vec_lagrange);

    VecDestroy(&vec_mu0);
    VecDestroy(&vec_mu1);
    VecDestroy(&vec_nest0);
    VecDestroy(&vec_nest1);
    return 0;
  }


180
181
182
183
184
185
  int petscMultMatFetiInterface(Mat mat, Vec x, Vec y)
  {
    FUNCNAME("petscMultMatFetiInterface()");

    double wtime = MPI::Wtime();

186
187
188
189
190
191
    Vec x_interface, x_lagrange, y_interface, y_lagrange;    
    VecNestGetSubVec(x, 0, &x_interface);
    VecNestGetSubVec(x, 1, &x_lagrange);
    VecNestGetSubVec(y, 0, &y_interface);
    VecNestGetSubVec(y, 1, &y_lagrange);

192
193
    void *ctx;
    MatShellGetContext(mat, &ctx);
194
    FetiData* data = static_cast<FetiData*>(ctx);
195
196


Thomas Witkowski's avatar
Thomas Witkowski committed
197
    // === Calculation of v_{B} ===
198

Thomas Witkowski's avatar
Thomas Witkowski committed
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
    // tmp_vec_b0 = J^{T} \lambda
    MatMultTranspose(*(data->mat_lagrange), x_lagrange, data->tmp_vec_b0);
    // tmp_vec_b1 = A_{B\Gamma} u_{\Gamma}
    MatMult(data->subSolver->getMatInteriorCoarse(1), x_interface, data->tmp_vec_b1);
    // tmp_vec_b0 = A_{B\Gamma} u_{\Gamma} + J^{T} \lambda
    VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1);


    // == Calculation of v_{\Pi}
    
    // tmp_vec_primal0 = A_{\Pi\Gamma} u_{\Gamma}
    MatMult(data->subSolver->getMatCoarse(0, 1), x_interface, data->tmp_vec_primal0);
    

    // === Calculate action of FETI-DP operator ===
214

Thomas Witkowski's avatar
Thomas Witkowski committed
215
216
217
218
    double wtime01 = MPI::Wtime();
    // tmp_vec_b0 = A_{BB}^{-1} v_{B}
    data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
219

Thomas Witkowski's avatar
Thomas Witkowski committed
220
    // tmp_vec_primal1 = A_{\Pi B} A_{BB}^{-1} v_{B}
221
    MatMult(data->subSolver->getMatCoarseInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
222
223
224
225
	    data->tmp_vec_b0, data->tmp_vec_primal1);

    // tmp_vec_primal0 = v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B}
    VecAXPY(data->tmp_vec_primal0, -1.0, data->tmp_vec_primal1);
226
227

    wtime01 = MPI::Wtime();
Thomas Witkowski's avatar
Thomas Witkowski committed
228
229
    // tmp_vec_primal0 = S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
    KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal0, data->tmp_vec_primal0);
230
231
    FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
232
    // tmp_vec_b1 = A_{B\Pi} S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
233
    MatMult(data->subSolver->getMatInteriorCoarse(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
234
	    data->tmp_vec_primal0, data->tmp_vec_b1);
235
236

    wtime01 = MPI::Wtime();
Thomas Witkowski's avatar
Thomas Witkowski committed
237
238
    // tmp_vec_b1 = A_{BB}^{-1} A_{B\Pi} S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
    data->subSolver->solveGlobal(data->tmp_vec_b1, data->tmp_vec_b1);
239
240
    FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);

Thomas Witkowski's avatar
Thomas Witkowski committed
241
242
243
244
245
246
247
248
249
250
251
252
253
254
    // tmp_vec_b0 = A_{BB}^{-1} v_{B} - A_{BB}^{-1} A_{B\Pi} S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
    VecAXPY(data->tmp_vec_b0, -1.0, data->tmp_vec_b1);


    // === Calculate projection to interface and constraint variable ===

    // y_interface = A_{\Gamma B} tmp_vec_b0
    MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, y_interface);

    // tmp_vec_primal0 = S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
    // tmp_vec_interface = A_{\Gamma \Pi} tmp_vec_primal0
    MatMult(data->subSolver->getMatCoarse(1, 0), data->tmp_vec_primal0, data->tmp_vec_interface);
    // y_interface = A_{\Gamma B} tmp_vec_b0 + A_{\Gamma \Pi} tmp_vec_primal0
    VecAXPY(y_interface, 1.0, data->tmp_vec_interface);
255

Thomas Witkowski's avatar
Thomas Witkowski committed
256
257
    // y_lagrange = J tmp_vec_b0
    MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y_lagrange);
258
259
260
261
262
263
264

    FetiTimings::fetiSolve += (MPI::Wtime() - wtime);

    return 0;
  }


265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
  PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
  {
    double wtime = MPI::Wtime();

    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);

    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);


    // === Restriction of the B nodes to the boundary nodes. ===

    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);

    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];

    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD - K_DI inv(K_II) K_ID ===

    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);

    MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
    KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
    MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);

    VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);


    // === Prolongation from local dual nodes to B nodes.

    VecGetArray(data->tmp_vec_b, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];

    VecRestoreArray(data->tmp_vec_b, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);

    FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);

    return 0;
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
330
  PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec xvec, Vec yvec)
331
332
333
334
335
336
  {
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);

Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
    // Multiply with scaled Lagrange constraint matrix.
    MatMultTranspose(*(data->mat_lagrange_scaled), xvec, data->tmp_vec_b0);


    // === Restriction of the B nodes to the boundary nodes. ===

    int nLocalB;
    int nLocalDuals;
    VecGetLocalSize(data->tmp_vec_b0, &nLocalB);
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);

    PetscScalar *local_b, *local_duals;
    VecGetArray(data->tmp_vec_b0, &local_b);
    VecGetArray(data->tmp_vec_duals0, &local_duals);

    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];

    VecRestoreArray(data->tmp_vec_b0, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // === K_DD ===

    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);


    // === Prolongation from local dual nodes to B nodes.

    VecGetArray(data->tmp_vec_b0, &local_b);
    VecGetArray(data->tmp_vec_duals1, &local_duals);

    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];

    VecRestoreArray(data->tmp_vec_b0, &local_b);
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);


    // Multiply with scaled Lagrange constraint matrix.
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b0, yvec);

    return 0;
  }


  PetscErrorCode petscApplyFetiInterfaceLumpedPrecon(PC pc, Vec xvec, Vec yvec)
  {
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
387
    FUNCNAME("precon");
Thomas Witkowski's avatar
Thomas Witkowski committed
388
389
390
391
392
393
394
395
396
397
398
399
    // Get data for the preconditioner
    void *ctx;
    PCShellGetContext(pc, &ctx);
    FetiInterfaceLumpedPreconData* data = 
      static_cast<FetiInterfaceLumpedPreconData*>(ctx);
    
    Vec x_interface, x_lagrange, y_interface, y_lagrange;    
    VecNestGetSubVec(xvec, 0, &x_interface);
    VecNestGetSubVec(xvec, 1, &x_lagrange);
    VecNestGetSubVec(yvec, 0, &y_interface);
    VecNestGetSubVec(yvec, 1, &y_lagrange);

400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
    bool useMassMatrix = false;
    Parameters::get("lp->mass matrix", useMassMatrix);
    if (useMassMatrix) {
      KSPSolve(data->ksp_mass, x_interface, y_interface);
    } else {
      VecCopy(x_interface, y_interface);
      double scaleFactor = 1.0;
      Parameters::get("lp->scal", scaleFactor);
      bool mult = false;
      Parameters::get("lp->mult", mult);
      if (mult)
	VecPointwiseMult(y_interface, x_interface, data->h2vec);
      if (scaleFactor != 0.0)
	VecScale(y_interface, scaleFactor);
    }
415

Thomas Witkowski's avatar
Thomas Witkowski committed
416
    MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0);   
417
418

    // === Restriction of the B nodes to the boundary nodes. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
419
    
420
421
    int nLocalB;
    int nLocalDuals;
Thomas Witkowski's avatar
Thomas Witkowski committed
422
    VecGetLocalSize(data->tmp_vec_b0, &nLocalB);
423
    VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
Thomas Witkowski's avatar
Thomas Witkowski committed
424
    
425
    PetscScalar *local_b, *local_duals;
Thomas Witkowski's avatar
Thomas Witkowski committed
426
    VecGetArray(data->tmp_vec_b0, &local_b);
427
    VecGetArray(data->tmp_vec_duals0, &local_duals);
Thomas Witkowski's avatar
Thomas Witkowski committed
428
    
429
430
431
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_duals[it->second] = local_b[it->first];
Thomas Witkowski's avatar
Thomas Witkowski committed
432
    
Thomas Witkowski's avatar
Thomas Witkowski committed
433
    VecRestoreArray(data->tmp_vec_b0, &local_b);
434
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
Thomas Witkowski's avatar
Thomas Witkowski committed
435
436
    
    
437
    // === K_DD ===
Thomas Witkowski's avatar
Thomas Witkowski committed
438
    
439
    MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
    
    
442
    // === Prolongation from local dual nodes to B nodes.
Thomas Witkowski's avatar
Thomas Witkowski committed
443
    
Thomas Witkowski's avatar
Thomas Witkowski committed
444
    VecGetArray(data->tmp_vec_b0, &local_b);
445
    VecGetArray(data->tmp_vec_duals1, &local_duals);
Thomas Witkowski's avatar
Thomas Witkowski committed
446
    
447
448
449
    for (map<int, int>::iterator it = data->localToDualMap.begin();
	 it != data->localToDualMap.end(); ++it)
      local_b[it->first] = local_duals[it->second];
Thomas Witkowski's avatar
Thomas Witkowski committed
450
    
Thomas Witkowski's avatar
Thomas Witkowski committed
451
    VecRestoreArray(data->tmp_vec_b0, &local_b);
452
    VecRestoreArray(data->tmp_vec_duals0, &local_duals);
Thomas Witkowski's avatar
Thomas Witkowski committed
453
        
454
    // Multiply with scaled Lagrange constraint matrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
455
    MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b0, y_lagrange);
Thomas Witkowski's avatar
Thomas Witkowski committed
456
  
457
458
459
    return 0;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
460

461
}