C++ LAPACK::GESV方法代码示例

MxDimMatrix<T, DIM> MxDimMatrix<T, DIM>::inv() const {

  MxDimMatrix<T, DIM> thisT = this->transpose(); 
  MxDimMatrix<T, DIM> solnT(MxDimVector<T, DIM>(1));

  int info;
  int ipiv[DIM];
  Teuchos::LAPACK<int, T> lapack;
  lapack.GESV(DIM, DIM, &thisT(0, 0), DIM, ipiv, &solnT(0, 0), DIM, &info);
  
  if (info != 0)
    std::cout << "MxDimMatrix::inv(): error in lapack routine. Return code: " << info << ".\n";

  return solnT.transpose();
}

//.........这里部分代码省略.........

              // compute neumann b.c. contributions and adjust rhs
              sideCub->getCubature(cub_points_side, cub_weights_side);
              for (unsigned i=0; i<numSides; i++) {
                // compute geometric cell information
                CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell, cub_points_side, sideDim, (int)i, cell);
                CellTools<double>::setJacobian(jacobian_side_refcell, cub_points_side_refcell, cell_nodes, cell);
                CellTools<double>::setJacobianDet(jacobian_det_side_refcell, jacobian_side_refcell);

                // compute weighted face measure
                FunctionSpaceTools::computeFaceMeasure<double>(weighted_measure_side_refcell,
                                                               jacobian_side_refcell,
                                                               cub_weights_side,
                                                               i,
                                                               cell);

                // tabulate values of basis functions at side cubature points, in the reference parent cell domain
                basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
                // transform 
                FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
                                                                value_of_basis_at_cub_points_side_refcell);

                // multiply with weighted measure
                FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
                                                            weighted_measure_side_refcell,
                                                            transformed_value_of_basis_at_cub_points_side_refcell);

                // compute Neumann data
                // map side cubature points in reference parent cell domain to physical space
                CellTools<double>::mapToPhysicalFrame(cub_points_side_physical, cub_points_side_refcell, cell_nodes, cell);
                // now compute data
                neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
                        cell, (int)i, x_order, y_order, z_order);

                FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
                                                      neumann_data_at_cub_points_side_physical,
                                                      weighted_transformed_value_of_basis_at_cub_points_side_refcell,
                                                      COMP_BLAS);

                // adjust RHS
                RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
              }
              ///////////////////////////////

              /////////////////////////////
              // Solution of linear system:
              int info = 0;
              Teuchos::LAPACK<int, double> solver;
              solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
              /////////////////////////////

              ////////////////////////
              // Building interpolant:
              // evaluate basis at interpolation points
              basis->getValues(value_of_basis_at_interp_points_ref, interp_points_ref, OPERATOR_VALUE);
              // transform values of basis functions 
              FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points_ref,
                                                              value_of_basis_at_interp_points_ref);
              FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points_ref);
              ////////////////////////

              /******************* END COMPUTATION ***********************/
          
              RealSpaceTools<double>::subtract(interpolant, exact_solution);

              *outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
                         << x_order << ", " << y_order << ", " << z_order
                         << ") and finite element interpolant of order " << basis_order << ": "
                         << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
                            RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";

              if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
                  RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
                *outStream << "\n\nPatch test failed for solution polynomial order ("
                           << x_order << ", " << y_order << ", " << z_order << ") and basis order " << basis_order << "\n\n";
                errorFlag++;
              }
            } // end for basis_order
          } // end for z_order
        } // end for y_order
      } // end for x_order
    } // end for ptype

  }
  // Catch unexpected errors
  catch (std::logic_error err) {
    *outStream << err.what() << "\n\n";
    errorFlag = -1000;
  };

  if (errorFlag != 0)
    std::cout << "End Result: TEST FAILED\n";
  else
    std::cout << "End Result: TEST PASSED\n";

  // reset format state of std::cout
  std::cout.copyfmt(oldFormatState);
 Kokkos::finalize();
  return errorFlag;
}

//.........这里部分代码省略.........
              
              // get normal direction, this has the edge weight factored into it already
              CellTools<double>::getReferenceSideNormal(side_normal , 
                                                        (int)side_cur,cell );

              // v.n at cub points on side
              vectorBasis->getValues(value_of_v_basis_at_cub_points_side ,
                                    cub_points_side_refcell ,
                                    OPERATOR_VALUE );

              for (int i=0;i<numVectorFields;i++) {
                for (int j=0;j<numCubPointsSide;j++) {
                  n_of_v_basis_at_cub_points_side(i,j) = 0.0;
                  for (int k=0;k<cellDim;k++) {
                    n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) * 
                      value_of_v_basis_at_cub_points_side(i,j,k);
                  }
                } 
              }
              
              cub_weights_side.resize(1,numCubPointsSide);
              FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side,
                                                          cub_weights_side,
                                                          n_of_v_basis_at_cub_points_side);
              cub_weights_side.resize(numCubPointsSide);
              
              FunctionSpaceTools::integrate<double>(rhs_vector_vec,
                                                    diri_data_at_cub_points_side,
                                                    w_n_of_v_basis_at_cub_points_side,
                                                    COMP_BLAS,
                                                    false);

              for (int i=0;i<numVectorFields;i++) {
                rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i);
              }
              
            }
            
            // solve linear system
            int info = 0;
            Teuchos::LAPACK<int, double> solver;
            solver.GESV(numTotalFields, 1, &fe_matrix(0,0,0), numTotalFields, &ipiv(0), &rhs_and_soln_vec(0,0), 
                        numTotalFields, &info);
            
            // compute interpolant; the scalar entries are last
            scalarBasis->getValues(value_of_s_basis_at_interp_points,
                                  interp_points_ref,
                                  OPERATOR_VALUE);
            for (int pt=0;pt<numInterpPoints;pt++) {
              interpolant(0,pt)=0.0;
              for (int i=0;i<numScalarFields;i++) {
                interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i)
                  * value_of_s_basis_at_interp_points(i,pt);
              }
            }
            
            interp_points_ref.resize(1,numInterpPoints,cellDim);
            // get exact solution for comparison
            FieldContainer<double> exact_solution(1,numInterpPoints);
            u_exact( exact_solution , interp_points_ref , x_order, y_order, z_order);
            interp_points_ref.resize(numInterpPoints,cellDim);

            RealSpaceTools<double>::add(interpolant,exact_solution);
            
            double nrm= RealSpaceTools<double>::vectorNorm(&interpolant(0,0),interpolant.dimension(1), NORM_TWO);

            *outStream << "\nNorm-2 error between scalar components of exact solution of order ("
                       << x_order << ", " << y_order << ", " << z_order
                       << ") and finite element interpolant of order " << basis_order << ": "
                       << nrm << "\n";

            if (nrm > zero) {
              *outStream << "\n\nPatch test failed for solution polynomial order ("
                         << x_order << ", " << y_order << ", " << z_order << ") and basis order (scalar, vector)  ("
                         << basis_order << ", " << basis_order+1 << ")\n\n";
              errorFlag++;
            }
            
          }
        }
      }
    }
    
  }
  
  catch (std::logic_error err) {
    *outStream << err.what() << "\n\n";
    errorFlag = -1000;
  };
  
  if (errorFlag != 0)
    std::cout << "End Result: TEST FAILED\n";
  else
    std::cout << "End Result: TEST PASSED\n";
  
  // reset format state of std::cout
  std::cout.copyfmt(oldFormatState);
  Kokkos::finalize();
  return errorFlag;
}

//.........这里部分代码省略.........
                                                rhs_at_cub_points_cell_physical,
                                                weighted_transformed_value_of_basis_at_cub_points_cell,
                                                COMP_BLAS);

          // compute neumann b.c. contributions and adjust rhs
          sideCub->getCubature(cub_points_side, cub_weights_side);
          for (unsigned i=0; i<numSides; i++) {
            // compute geometric cell information
            CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell, cub_points_side, sideDim, (int)i, cell);
            CellTools<double>::setJacobian(jacobian_side_refcell, cub_points_side_refcell, cell_nodes[pcell], cell);
            CellTools<double>::setJacobianDet(jacobian_det_side_refcell, jacobian_side_refcell);

            // compute weighted edge measure
            FunctionSpaceTools::computeEdgeMeasure<double>(weighted_measure_side_refcell,
                                                           jacobian_side_refcell,
                                                           cub_weights_side,
                                                           i,
                                                           cell);

            // tabulate values of basis functions at side cubature points, in the reference parent cell domain
            basis->getValues(value_of_basis_at_cub_points_side_refcell, cub_points_side_refcell, OPERATOR_VALUE);
            // transform
            FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points_side_refcell,
                                                            value_of_basis_at_cub_points_side_refcell);

            // multiply with weighted measure
            FunctionSpaceTools::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_side_refcell,
                                                        weighted_measure_side_refcell,
                                                        transformed_value_of_basis_at_cub_points_side_refcell);

            // compute Neumann data
            // map side cubature points in reference parent cell domain to physical space
            CellTools<double>::mapToPhysicalFrame(cub_points_side_physical, cub_points_side_refcell, cell_nodes[pcell], cell);
            // now compute data
            neumann(neumann_data_at_cub_points_side_physical, cub_points_side_physical, jacobian_side_refcell,
                    cell, (int)i, x_order, y_order);

            FunctionSpaceTools::integrate<double>(neumann_fields_per_side,
                                                  neumann_data_at_cub_points_side_physical,
                                                  weighted_transformed_value_of_basis_at_cub_points_side_refcell,
                                                  COMP_BLAS);

            // adjust RHS
            RealSpaceTools<double>::add(rhs_and_soln_vector, neumann_fields_per_side);;
          }
          ///////////////////////////////

          /////////////////////////////
          // Solution of linear system:
          int info = 0;
          Teuchos::LAPACK<int, double> solver;
          solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vector[0], numFields, &info);
          /////////////////////////////

          ////////////////////////
          // Building interpolant:
          // evaluate basis at interpolation points
          basis->getValues(value_of_basis_at_interp_points, interp_points_ref, OPERATOR_VALUE);
          // transform values of basis functions
          FunctionSpaceTools::HGRADtransformVALUE<double>(transformed_value_of_basis_at_interp_points,
                                                          value_of_basis_at_interp_points);
          FunctionSpaceTools::evaluate<double>(interpolant, rhs_and_soln_vector, transformed_value_of_basis_at_interp_points);
          ////////////////////////

          /******************* END COMPUTATION ***********************/
      
          RealSpaceTools<double>::subtract(interpolant, exact_solution);

          *outStream << "\nRelative norm-2 error between exact solution polynomial of order ("
                     << x_order << ", " << y_order << ") and finite element interpolant of order " << basis_order << ": "
                     << RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
                        RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) << "\n";

          if (RealSpaceTools<double>::vectorNorm(&interpolant[0], interpolant.dimension(1), NORM_TWO) /
              RealSpaceTools<double>::vectorNorm(&exact_solution[0], exact_solution.dimension(1), NORM_TWO) > zero) {
            *outStream << "\n\nPatch test failed for solution polynomial order ("
                       << x_order << ", " << y_order << ") and basis order " << basis_order << "\n\n";
            errorFlag++;
          }
        } // end for y_order
      } // end for x_order
    } // end for pcell

  }
  // Catch unexpected errors
  catch (std::logic_error err) {
    *outStream << err.what() << "\n\n";
    errorFlag = -1000;
  };

  if (errorFlag != 0)
    std::cout << "End Result: TEST FAILED\n";
  else
    std::cout << "End Result: TEST PASSED\n";

  // reset format state of std::cout
  std::cout.copyfmt(oldFormatState);

  return errorFlag;
}

//.........这里部分代码省略.........
                    if (verbose && p < 5) {
                      *outStream << " - Element matrix and rhs, iel = " << iel << "\n";
                      *outStream << std::showpos;
                      for (int i=0;i<nbf;++i) {
                        for (int j=0;j<nbf;++j)
                          *outStream << MatVal(A, i, j) << " ";
                        *outStream << ":: " << MatVal(b, i, 0) << "\n";
                      }
                      *outStream << std::noshowpos;
                    }

                    // Step 1.2: assemble high order elements
                    assemble_element_matrix_and_rhs(A_asm, b_asm,
                                                    A, b,
                                                    local2global[iel],
                                                    nnodes_per_element);
                  }

                  if (verbose && p < 5) {
                    *outStream << " - Assembled element matrix and rhs -\n";
                    *outStream << std::showpos;
                    for (int i=0;i<ndofs;++i) {
                      for (int j=0;j<ndofs;++j)
                        *outStream << MatVal(A_asm, i, j) << " ";
                      *outStream << ":: " << MatVal(b_asm, i, 0) << "\n";
                    }
                    *outStream << std::noshowpos;
                  }

                  // Step 2: solve the system of equations
                  int info = 0;
                  Teuchos::LAPACK<int,value_type> lapack;
                  FieldContainer<int> ipiv(ndofs);
                  lapack.GESV(ndofs, 1, &A_asm(0,0,0), ndofs, &ipiv(0,0), &b_asm(0,0), ndofs, &info);
                  TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::runtime_error,
                                              ">>> ERROR (Intrepid::HGRAD_TRI_Cn::Test 04): " \
                                              "LAPACK solve fails");

                  // Step 3: construct interpolant and check solutions
                  magnitude_type interpolation_error = 0, solution_norm =0;
                  for (int iel=0;iel<nelement;++iel) {
                    retrieve_element_solution(b,
                                              b_asm,
                                              local2global[iel],
                                              nnodes_per_element);

                    if (verbose && p < 5) {
                      *outStream << " - Element solution, iel = " << iel << "\n";
                      *outStream << std::showpos;
                      for (int i=0;i<nbf;++i) {
                        *outStream << MatVal(b, i, 0) << "\n";
                      }
                      *outStream << std::noshowpos;
                    }

                    magnitude_type
                      element_interpolation_error = 0,
                      element_solution_norm = 0;

                    Orientation ort = Orientation::getOrientation(cell, element[iel]);

                    // set element nodal coordinates
                    fill_cell_nodes(cell_nodes, 
                                    nodes, 
                                    element[iel], 
                                    nvert, ndim);

//.........这里部分代码省略.........
  Teuchos::RCP<NOX::Abstract::MultiVector> tmp = 
    result_null.clone(NOX::ShapeCopy);
  status = group->computeDwtJnDxMulti(result_null, *nullVector, *tmp);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							   finalStatus,
							   callingFunction);

  // compute [F 0 0] - (Jv)_x^T[A B u]
  tmp->update(1.0, input_x, -1.0);

  // verify underlying Jacobian is valid
  if (!group->isJacobian()) {
    status = group->computeJacobian();
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // Solve  |J^T v||C D E| = |F - (Jv)_x^T A  -(Jv)_x^T B  -(Jv)_x^T u|
  //        |u^T 0||c d e|   |         0             0            0   |
  status = 
    transposeBorderedSolver->applyInverseTranspose(params, 
						   tmp.get(), 
						   NULL, 
						   result_x,
						   tmp_mat_2);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
							   callingFunction);
  Teuchos::RCP<NOX::Abstract::MultiVector> C = 
    result_x.subView(index_input);
  Teuchos::RCP<NOX::Abstract::MultiVector> D = 
    result_x.subView(index_dp);
  Teuchos::RCP<NOX::Abstract::MultiVector> E = 
    result_x.subView(index_null);
  double d = tmp_mat_2(0, m);
  double e = tmp_mat_2(0, m+1);

  // compute (Jv)_p^T*[A B u]
  NOX::Abstract::MultiVector::DenseMatrix t1(1,m+2);
  result_null.multiply(1.0, *dJndp, t1);

  // compute f_p^T*[C D E]
  NOX::Abstract::MultiVector::DenseMatrix t2(1,m+2);
  result_x.multiply(1.0, *dfdp, t2);

  // compute f_p^T*u
  double fptu = uVector->innerProduct((*dfdp)[0]);

  // Fill coefficient arrays
  double M[9];
  M[0] = st;   M[1] =  -e;   M[2] = t1(0,m+1) + t2(0,m+1);
  M[3] = 0.0;  M[4] =   st;  M[5] = fptu;
  M[6] = -b;   M[7] =  -d;   M[8] = t1(0,m) + t2(0,m);

  // Compute RHS
  double *R = new double[3*m];
  for (int i=0; i<m; i++) {
    R[3*i]   = tmp_mat_1(0,i);
    R[3*i+1] = tmp_mat_2(0,i);
    R[3*i+2] = result_param(0,i) - t1(0,i) - t2(0,i);
  }

  // Solve M*P = R
  int three = 3;
  int piv[3];
  int info;
  Teuchos::LAPACK<int,double> L;
  L.GESV(three, m, M, three, piv, R, three, &info);
  if (info != 0) {
    globalData->locaErrorCheck->throwError(
				    callingFunction,
				    "Solve of 3x3 coefficient matrix failed!");
    return NOX::Abstract::Group::Failed;
  }

  NOX::Abstract::MultiVector::DenseMatrix alpha(1,m);
  NOX::Abstract::MultiVector::DenseMatrix beta(1,m);
  for (int i=0; i<m; i++) {
    alpha(0,i)        = R[3*i];
    beta(0,i)         = R[3*i+1];
    result_param(0,i) = R[3*i+2];
  }

  // compute A = A + B*z + alpha*u (remember A is a sub-view of result_null)
  A->update(Teuchos::NO_TRANS, 1.0, *B, result_param, 1.0);
  A->update(Teuchos::NO_TRANS, 1.0, *uMultiVector, alpha, 1.0);

  // compute C = C + D*z + alpha*E + beta*u 
  // (remember C is a sub-view of result_x)
  C->update(Teuchos::NO_TRANS, 1.0, *D, result_param, 1.0);
  C->update(Teuchos::NO_TRANS, 1.0, *E, alpha, 1.0);
  C->update(Teuchos::NO_TRANS, 1.0, *uMultiVector, beta, 1.0);

  delete [] R;

  return finalStatus;
}

//.........这里部分代码省略.........
  // verify underlying Jacobian is valid
  if (!group->isJacobian()) {
    status = group->computeJacobian();
    finalStatus = 
      globalData->locaErrorCheck->combineAndCheckReturnTypes(status, 
							     finalStatus,
							     callingFunction);
  }

  // Solve  |J   u||D E K L| = |G-(Jv)_xA  d(Jv)/dp-(Jv)_xB  -(Jv)_xC -(Jv)_xv|
  //        |v^T 0||d e k l|   |    0             0              0        0   |
  status = borderedSolver->applyInverse(params, tmp.get(), NULL, result_null,
					tmp_mat_2);
  finalStatus = 
    globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
							   callingFunction);
  Teuchos::RCP<NOX::Abstract::MultiVector> D = 
    result_null.subView(index_input);
  Teuchos::RCP<NOX::Abstract::MultiVector> E = 
    result_null.subView(index_dp);
  Teuchos::RCP<NOX::Abstract::MultiVector> K = 
    result_null.subView(index_s);
  Teuchos::RCP<NOX::Abstract::MultiVector> L = 
    result_null.subView(index_null);
  double e = tmp_mat_2(0, m);
  double k = tmp_mat_2(0, m+1);
  double l = tmp_mat_2(0, m+2);

  double ltE = pfGroup->lTransNorm((*E)[0]);
  double ltK = pfGroup->lTransNorm((*K)[0]);
  double ltL = pfGroup->lTransNorm((*L)[0]);
  double ltv = pfGroup->lTransNorm(*nullVector);
  double ipv = group->innerProduct(*nullVector, *asymVector);
  double ipB = group->innerProduct((*B)[0], *asymVector);
  double ipC = group->innerProduct((*C)[0], *asymVector);

  // Fill coefficient arrays
  double M[16];
  M[0]  = sigma; M[1]  = -l;     M[2]  =  ipv; M[3]  =  ltL;
  M[4]  = 0.0;   M[5]  =  sigma; M[6]  =  0.0; M[7]  =  ltv;
  M[8]  = b;     M[9]  =  e;     M[10] = -ipB; M[11] = -ltE;
  M[12] = c;     M[13] =  k;     M[14] = -ipC; M[15] = -ltK;

  // compute s - <A,psi>
  NOX::Abstract::MultiVector::DenseMatrix tmp_mat_3(1, m);
  group->innerProduct(*asymMultiVector, *A, tmp_mat_3);
  tmp_mat_3 -= input_slack;
  tmp_mat_3.scale(-1.0);

  // compute h - phi^T D
  NOX::Abstract::MultiVector::DenseMatrix tmp_mat_4(1, m);
  pfGroup->lTransNorm(*D, tmp_mat_4);
  tmp_mat_4 -= input_param;
  tmp_mat_4.scale(-1.0);

  double *R = new double[4*m];
  for (int i=0; i<m; i++) {
    R[4*i]   = tmp_mat_1(0,i);
    R[4*i+1] = tmp_mat_2(0,i);
    R[4*i+2] = tmp_mat_3(0,i);
    R[4*i+3] = tmp_mat_4(0,i);
  }

  // Solve M*P = R
  int piv[4];
  int info;
  Teuchos::LAPACK<int,double> dlapack;
  dlapack.GESV(4, m, M, 4, piv, R, 4, &info);
  if (info != 0) {
    globalData->locaErrorCheck->throwError(
				    callingFunction,
				    "Solve of 4x4 coefficient matrix failed!");
    return NOX::Abstract::Group::Failed;
  }

  NOX::Abstract::MultiVector::DenseMatrix alpha(1,m);
  NOX::Abstract::MultiVector::DenseMatrix beta(1,m);
  for (int i=0; i<m; i++) {
    alpha(0,i)        = R[4*i];
    beta(0,i)         = R[4*i+1];
    result_param(0,i) = R[4*i+2];
    result_slack(0,i) = R[4*i+3];
  }

  // compute A = A - B*z -C*w + v*alpha (remember A is a sub-view of result_x)
  A->update(Teuchos::NO_TRANS, -1.0, *B, result_param, 1.0);
  A->update(Teuchos::NO_TRANS, -1.0, *C, result_slack, 1.0);
  A->update(Teuchos::NO_TRANS, 1.0, *nullMultiVector, alpha, 1.0);

  // compute D = D - E*z - K*w + L*alpha + v*beta 
  // (remember D is a sub-view of result_null)
  D->update(Teuchos::NO_TRANS, -1.0, *E, result_param, 1.0);
  D->update(Teuchos::NO_TRANS, -1.0, *K, result_slack, 1.0);
  D->update(Teuchos::NO_TRANS, 1.0, *L, alpha, 1.0);
  D->update(Teuchos::NO_TRANS, 1.0, *nullMultiVector, beta, 1.0);

  delete [] R;

  return finalStatus;
}