List of usage examples for org.apache.commons.math.linear FieldMatrix setEntry
void setEntry(int row, int column, T value) throws MatrixIndexException;
From source file:rb.app.ComplexRBSystem.java
@Override public double RB_solve(int N) { current_N = N;/*from w ww. jav a 2 s. co m*/ theta_a = complex_eval_theta_q_a(); if (N > get_n_basis_functions()) { throw new RuntimeException( "ERROR: N cannot be larger than the number " + "of basis functions in RB_solve"); } if (N == 0) { throw new RuntimeException("ERROR: N must be greater than 0 in RB_solve"); } // Assemble the RB system FieldMatrix<Complex> RB_system_matrix_N = new Array2DRowFieldMatrix<Complex>((new Complex(0, 0)).getField(), N, N); for (int q_a = 0; q_a < get_Q_a(); q_a++) { RB_system_matrix_N = RB_system_matrix_N .add(RB_A_q_vector[q_a].getSubMatrix(0, N - 1, 0, N - 1).scalarMultiply(theta_a.getEntry(q_a))); // scalarMultiply(complex_eval_theta_q_a(q_a) ) ); } // Assemble the RB rhs FieldVector<Complex> RB_rhs_N = new ArrayFieldVector<Complex>(N, new Complex(0d, 0d)); for (int q_f = 0; q_f < get_Q_f(); q_f++) { // Note getSubVector takes an initial index and the number of entries // i.e. the interface is a bit different to getSubMatrix RB_rhs_N = RB_rhs_N.add(RB_F_q_vector[q_f].getSubVector(0, N).mapMultiply(complex_eval_theta_q_f(q_f))); } // Solve the linear system by Gaussian elimination RB_solution = new ArrayFieldVector<Complex>(N, new Complex(0., 0.)); for (int j = 1; j < N; j++) for (int i = j; i < N; i++) { Complex m = RB_system_matrix_N.getEntry(i, j - 1).divide(RB_system_matrix_N.getEntry(j - 1, j - 1)); for (int k = 0; k < N; k++) RB_system_matrix_N.setEntry(i, k, RB_system_matrix_N.getEntry(i, k) .subtract(RB_system_matrix_N.getEntry(j - 1, k).multiply(m))); RB_rhs_N.setEntry(i, RB_rhs_N.getEntry(i).subtract(m.multiply(RB_rhs_N.getEntry(j - 1)))); } RB_solution.setEntry(N - 1, RB_rhs_N.getEntry(N - 1).divide(RB_system_matrix_N.getEntry(N - 1, N - 1))); for (int j = N - 2; j >= 0; j--) { Complex m = new Complex(0., 0.); for (int i = j + 1; i < N; i++) m = m.add(RB_system_matrix_N.getEntry(j, i).multiply(RB_solution.getEntry(i))); RB_solution.setEntry(j, (RB_rhs_N.getEntry(j).subtract(m)).divide(RB_system_matrix_N.getEntry(j, j))); } // Evaluate the dual norm of the residual for RB_solution_vector double epsilon_N = compute_residual_dual_norm(N); // Get lower bound for coercivity constant double alpha_LB = get_SCM_lower_bound(); // If SCM lower bound is negative if (alpha_LB < 0) { // Get an upper bound instead alpha_LB = get_SCM_upper_bound(); } // Store (absolute) error bound double abs_error_bound = epsilon_N / residual_scaling_denom(alpha_LB); // Compute the norm of RB_solution /* RealMatrix RB_inner_product_matrix_N = RB_inner_product_matrix.getSubMatrix(0, N-1, 0, N-1); */ double RB_solution_norm = 0.0d; for (int i = 0; i < N; i++) RB_solution_norm += ((RB_solution.getEntry(i)).multiply((RB_solution.getEntry(i)).conjugate())) .getReal(); RB_solution_norm = Math.sqrt(RB_solution_norm); // Now compute the outputs and associated errors FieldVector<Complex> RB_output_vector_N = new ArrayFieldVector<Complex>(N, new Complex(0d, 0d)); for (int i = 0; i < get_n_outputs(); i++) { RB_outputs[i] = new Complex(0., 0.); RB_output_vector_N = (RB_output_vectors[i][0].getSubVector(0, N)) .mapMultiply(complex_eval_theta_q_l(i, 0)); for (int q_l = 1; q_l < get_Q_l(i); q_l++) RB_output_vector_N = RB_output_vector_N.add( (RB_output_vectors[i][q_l].getSubVector(0, N)).mapMultiply(complex_eval_theta_q_l(i, q_l))); for (int j = 0; j < N; j++) RB_outputs[i] = RB_outputs[i] .add((RB_output_vector_N.getEntry(j).conjugate()).multiply((RB_solution.getEntry(j)))); RB_output_error_bounds[i] = new Complex(compute_output_dual_norm(i) * abs_error_bound, compute_output_dual_norm(i) * abs_error_bound); } cal_derived_output(); return (return_rel_error_bound ? abs_error_bound / RB_solution_norm : abs_error_bound); }
From source file:rb.ComplexRBSystem.java
@Override public double RB_solve(int N) { current_N = N;/* ww w. jav a2s. c o m*/ theta_a = complex_eval_theta_q_a(); if (N > getNBF()) { throw new RuntimeException( "ERROR: N cannot be larger than the number " + "of basis functions in RB_solve"); } if (N == 0) { throw new RuntimeException("ERROR: N must be greater than 0 in RB_solve"); } // Assemble the RB system FieldMatrix<Complex> RB_system_matrix_N = new Array2DRowFieldMatrix<Complex>((new Complex(0, 0)).getField(), N, N); for (int q_a = 0; q_a < getQa(); q_a++) { RB_system_matrix_N = RB_system_matrix_N .add(RB_A_q_vector[q_a].getSubMatrix(0, N - 1, 0, N - 1).scalarMultiply(theta_a.getEntry(q_a))); // scalarMultiply(complex_eval_theta_q_a(q_a) ) ); } // Assemble the RB rhs FieldVector<Complex> RB_rhs_N = new ArrayFieldVector<Complex>(N, new Complex(0d, 0d)); for (int q_f = 0; q_f < getQf(); q_f++) { // Note getSubVector takes an initial index and the number of // entries // i.e. the interface is a bit different to getSubMatrix RB_rhs_N = RB_rhs_N.add(RB_F_q_vector[q_f].getSubVector(0, N).mapMultiply(complex_eval_theta_q_f(q_f))); } // Solve the linear system by Gaussian elimination RB_solution = new ArrayFieldVector<Complex>(N, new Complex(0., 0.)); for (int j = 1; j < N; j++) for (int i = j; i < N; i++) { Complex m = RB_system_matrix_N.getEntry(i, j - 1).divide(RB_system_matrix_N.getEntry(j - 1, j - 1)); for (int k = 0; k < N; k++) RB_system_matrix_N.setEntry(i, k, RB_system_matrix_N.getEntry(i, k) .subtract(RB_system_matrix_N.getEntry(j - 1, k).multiply(m))); RB_rhs_N.setEntry(i, RB_rhs_N.getEntry(i).subtract(m.multiply(RB_rhs_N.getEntry(j - 1)))); } RB_solution.setEntry(N - 1, RB_rhs_N.getEntry(N - 1).divide(RB_system_matrix_N.getEntry(N - 1, N - 1))); for (int j = N - 2; j >= 0; j--) { Complex m = new Complex(0., 0.); for (int i = j + 1; i < N; i++) m = m.add(RB_system_matrix_N.getEntry(j, i).multiply(RB_solution.getEntry(i))); RB_solution.setEntry(j, (RB_rhs_N.getEntry(j).subtract(m)).divide(RB_system_matrix_N.getEntry(j, j))); } // Evaluate the dual norm of the residual for RB_solution_vector double epsilon_N = compute_residual_dual_norm(N); // Get lower bound for coercivity constant double alpha_LB = get_SCM_lower_bound(); // If SCM lower bound is negative if (alpha_LB < 0) { // Get an upper bound instead alpha_LB = get_SCM_upper_bound(); } // Store (absolute) error bound double abs_error_bound = epsilon_N / residual_scaling_denom(alpha_LB); // Compute the norm of RB_solution /* * RealMatrix RB_inner_product_matrix_N = * RB_inner_product_matrix.getSubMatrix(0, N-1, 0, N-1); */ double RB_solution_norm = 0.0d; for (int i = 0; i < N; i++) RB_solution_norm += ((RB_solution.getEntry(i)).multiply((RB_solution.getEntry(i)).conjugate())) .getReal(); RB_solution_norm = Math.sqrt(RB_solution_norm); // Now compute the outputs and associated errors FieldVector<Complex> RB_output_vector_N = new ArrayFieldVector<Complex>(N, new Complex(0d, 0d)); for (int i = 0; i < getNumOutputs(); i++) { RB_outputs[i] = new Complex(0., 0.); RB_output_vector_N = (RB_output_vectors[i][0].getSubVector(0, N)) .mapMultiply(complex_eval_theta_q_l(i, 0)); for (int q_l = 1; q_l < getQl(i); q_l++) RB_output_vector_N = RB_output_vector_N.add( (RB_output_vectors[i][q_l].getSubVector(0, N)).mapMultiply(complex_eval_theta_q_l(i, q_l))); for (int j = 0; j < N; j++) RB_outputs[i] = RB_outputs[i] .add((RB_output_vector_N.getEntry(j).conjugate()).multiply((RB_solution.getEntry(j)))); RB_output_error_bounds[i] = new Complex(compute_output_dual_norm(i, 0) // Zero // means // no // time // used // here * abs_error_bound, compute_output_dual_norm(i, 0) // Zero // means // no // time // used // here * abs_error_bound); } cal_derived_output(); return (return_rel_error_bound ? abs_error_bound / RB_solution_norm : abs_error_bound); }