List of usage examples for org.apache.commons.math.linear Array2DRowFieldMatrix Array2DRowFieldMatrix
public Array2DRowFieldMatrix(final Field<T> field, final int rowDimension, final int columnDimension) throws IllegalArgumentException
From source file:rb.app.ComplexRBSystem.java
@Override public void read_offline_data(String directory_name, boolean isAssetFile) throws IOException { isReal = false;/* w w w . j av a 2 s. c o m*/ HttpClient client = new DefaultHttpClient(); int buffer_size = 8192; // Find out dimension of the RB space { InputStreamReader isr; String dataString = directory_name + "/n_bfs.dat"; if (!isAssetFile) { // Read from server HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); isr = new InputStreamReader(response.getEntity().getContent()); } else { // Read from assets isr = new InputStreamReader(context.getAssets().open(dataString)); } BufferedReader reader = new BufferedReader(isr, buffer_size); String line = reader.readLine(); set_n_basis_functions(Integer.parseInt(line)); reader.close(); } Log.d(DEBUG_TAG, "Finished reading n_bfs.dat"); // Read in output data if (get_n_outputs() > 0) { // Get output dual norms { RB_output_vectors = (FieldVector<Complex>[][]) new FieldVector<?>[get_n_outputs()][]; output_dual_norms = new Complex[get_n_outputs()][]; for (int i = 0; i < get_n_outputs(); i++) { InputStreamReader isr1; String dataString1 = directory_name + "/output_" + String.format("%03d", i) + "_dual_norms.dat"; if (!isAssetFile) { // Read from server HttpGet request1 = new HttpGet(dataString1); HttpResponse response1 = client.execute(request1); isr1 = new InputStreamReader(response1.getEntity().getContent()); } else { // Read from assets isr1 = new InputStreamReader(context.getAssets().open(dataString1)); } BufferedReader reader1 = new BufferedReader(isr1, buffer_size); String line1 = reader1.readLine(); String[] dual_norms_tokens = line1.split(" "); int Q_l_hat = get_Q_l(i) * (get_Q_l(i) + 1) / 2; output_dual_norms[i] = new Complex[Q_l_hat]; for (int q = 0; q < Q_l_hat; q++) { output_dual_norms[i][q] = new Complex(Double.parseDouble(dual_norms_tokens[q]), Double.parseDouble(dual_norms_tokens[Q_l_hat + q])); } RB_output_vectors[i] = (FieldVector<Complex>[]) new FieldVector<?>[get_Q_l(i)]; for (int q_l = 0; q_l < get_Q_l(i); q_l++) { // Now read in the RB output vectors InputStreamReader isr_i; String dataString_i = directory_name + "/output_" + String.format("%03d", i) + "_" + String.format("%03d", q_l) + ".dat"; if (!isAssetFile) { // Read from server HttpGet request_i = new HttpGet(dataString_i); HttpResponse response_i = client.execute(request_i); isr_i = new InputStreamReader(response_i.getEntity().getContent()); } else { // Read from assets isr_i = new InputStreamReader(context.getAssets().open(dataString_i)); } BufferedReader reader_i = new BufferedReader(isr_i, buffer_size); String line_i = reader_i.readLine(); String[] output_i_tokens = line_i.split(" "); RB_output_vectors[i][q_l] = new ArrayFieldVector<Complex>(get_n_basis_functions(), new Complex(0d, 0d)); for (int j = 0; j < get_n_basis_functions(); j++) { RB_output_vectors[i][q_l].setEntry(j, new Complex(Double.parseDouble(output_i_tokens[j]), Double.parseDouble(output_i_tokens[get_n_basis_functions() + j]))); } reader_i.close(); } reader1.close(); } } } Log.d(DEBUG_TAG, "Finished reading output data"); /* // Read in the inner product matrix { InputStreamReader isr; String dataString = directory_name + "/RB_inner_product_matrix.dat"; if(!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); isr = new InputStreamReader(response.getEntity() .getContent()); } else { // Read from assets isr = new InputStreamReader( context.getAssets().open(dataString)); } BufferedReader reader = new BufferedReader(isr,buffer_size); String line = reader.readLine(); String[] tokens = line.split(" "); // Set the size of the inner product matrix RB_inner_product_matrix = new Array2DRowRealMatrix(n_bfs, n_bfs); // Fill the matrix int count = 0; for (int i = 0; i < n_bfs; i++) for (int j = 0; j < n_bfs; j++) { RB_inner_product_matrix.setEntry(i, j, Double .parseDouble(tokens[count])); count++; } reader.close(); } Log.d(DEBUG_TAG, "Finished reading RB_inner_product_matrix.dat"); */ // Read in the F_q vectors RB_F_q_vector = (FieldVector<Complex>[]) new FieldVector<?>[get_Q_f()]; for (int q_f = 0; q_f < get_Q_f(); q_f++) { InputStreamReader isr; String dataString = directory_name + "/RB_F_" + String.format("%03d", q_f) + ".dat"; if (!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); isr = new InputStreamReader(response.getEntity().getContent()); } else { // Read from assets isr = new InputStreamReader(context.getAssets().open(dataString)); } BufferedReader reader = new BufferedReader(isr, buffer_size); String line = reader.readLine(); String[] tokens = line.split(" "); // Set the size of the inner product matrix RB_F_q_vector[q_f] = new ArrayFieldVector<Complex>(get_n_basis_functions(), new Complex(0d, 0d)); // Fill the vector for (int i = 0; i < get_n_basis_functions(); i++) { RB_F_q_vector[q_f].setEntry(i, new Complex(Double.parseDouble(tokens[i]), Double.parseDouble(tokens[get_n_basis_functions() + i]))); } reader.close(); } Log.d(DEBUG_TAG, "Finished reading RB_F_q data"); // Read in the A_q matrices RB_A_q_vector = (FieldMatrix<Complex>[]) new FieldMatrix<?>[get_Q_a()]; for (int q_a = 0; q_a < get_Q_a(); q_a++) { InputStreamReader isr; String dataString = directory_name + "/RB_A_" + String.format("%03d", q_a) + ".dat"; if (!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); isr = new InputStreamReader(response.getEntity().getContent()); } else { // Read from assets isr = new InputStreamReader(context.getAssets().open(dataString)); } BufferedReader reader = new BufferedReader(isr, buffer_size); String line = reader.readLine(); String[] tokens = line.split(" "); // Set the size of the inner product matrix RB_A_q_vector[q_a] = new Array2DRowFieldMatrix<Complex>((new Complex(0, 0)).getField(), get_n_basis_functions(), get_n_basis_functions()); // Fill the vector int count = 0; for (int i = 0; i < get_n_basis_functions(); i++) for (int j = 0; j < get_n_basis_functions(); j++) { RB_A_q_vector[q_a].setEntry(i, j, new Complex(Double.parseDouble(tokens[count]), Double.parseDouble(tokens[count + get_n_basis_functions() * get_n_basis_functions()]))); count++; } reader.close(); } Log.d(DEBUG_TAG, "Finished reading RB_A_q data"); // Read in F_q representor norm data { InputStreamReader isr; String dataString = directory_name + "/Fq_norms.dat"; if (!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); isr = new InputStreamReader(response.getEntity().getContent()); } else { // Read from assets isr = new InputStreamReader(context.getAssets().open(dataString)); } BufferedReader reader = new BufferedReader(isr, buffer_size); String line = reader.readLine(); String[] tokens = line.split(" "); // Declare the array int Q_f_hat = get_Q_f() * (get_Q_f() + 1) / 2; Fq_representor_norms = new Complex[Q_f_hat]; // Fill it for (int i = 0; i < Q_f_hat; i++) { Fq_representor_norms[i] = new Complex(Double.parseDouble(tokens[i * 2 + 0]), Double.parseDouble(tokens[i * 2 + 1])); } reader.close(); } Log.d(DEBUG_TAG, "Finished reading Fq_norms.dat"); // Read in Fq_Aq representor norm data { InputStreamReader isr; String dataString = directory_name + "/Fq_Aq_norms.dat"; if (!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); isr = new InputStreamReader(response.getEntity().getContent()); } else { // Read from assets isr = new InputStreamReader(context.getAssets().open(dataString)); } BufferedReader reader = new BufferedReader(isr, buffer_size); String line = reader.readLine(); String[] tokens = line.split(" "); // Declare the array Fq_Aq_representor_norms = new Complex[get_Q_f()][get_Q_a()][get_n_basis_functions()]; double[][][] Rdata = new double[get_Q_f()][get_Q_a()][get_n_basis_functions()]; double[][][] Idata = new double[get_Q_f()][get_Q_a()][get_n_basis_functions()]; // Fill it int count = 0; for (int q_f = 0; q_f < get_Q_f(); q_f++) for (int q_a = 0; q_a < get_Q_a(); q_a++) { for (int i = 0; i < get_n_basis_functions(); i++) { Rdata[q_f][q_a][i] = Double.parseDouble(tokens[count]); count++; } for (int i = 0; i < get_n_basis_functions(); i++) { Idata[q_f][q_a][i] = Double.parseDouble(tokens[count]); count++; } } for (int q_f = 0; q_f < get_Q_f(); q_f++) for (int q_a = 0; q_a < get_Q_a(); q_a++) for (int i = 0; i < get_n_basis_functions(); i++) Fq_Aq_representor_norms[q_f][q_a][i] = new Complex(Rdata[q_f][q_a][i], Idata[q_f][q_a][i]); reader.close(); } Log.d(DEBUG_TAG, "Finished reading Fq_Aq_norms.dat"); // Read in Aq_Aq representor norm data { // Declare the array int Q_a_hat = get_Q_a() * (get_Q_a() + 1) / 2; Aq_Aq_representor_norms = new Complex[Q_a_hat][get_n_basis_functions()][get_n_basis_functions()]; int count = 0; for (int i = 0; i < get_Q_a(); i++) for (int j = i; j < get_Q_a(); j++) { BinaryReader dis; String dataString = directory_name + "/Aq_Aq_" + String.format("%03d", i) + "_" + String.format("%03d", j) + "_norms.bin"; if (!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); dis = new BinaryReader(response.getEntity().getContent()); } else { // Read from assets dis = new BinaryReader(context.getAssets().open(dataString)); } double[][] Rdata = new double[get_n_basis_functions()][get_n_basis_functions()]; double[][] Idata = new double[get_n_basis_functions()][get_n_basis_functions()]; for (int k = 0; k < get_n_basis_functions(); k++) for (int l = 0; l < get_n_basis_functions(); l++) Rdata[k][l] = dis.ReadDouble(); for (int k = 0; k < get_n_basis_functions(); k++) for (int l = 0; l < get_n_basis_functions(); l++) Idata[k][l] = dis.ReadDouble(); for (int k = 0; k < get_n_basis_functions(); k++) for (int l = 0; l < get_n_basis_functions(); l++) Aq_Aq_representor_norms[count][k][l] = new Complex(Rdata[k][l], Idata[k][l]); count++; dis.close(); } } Log.d(DEBUG_TAG, "Finished reading Aq_Aq_norms.dat"); // Read calN number if (get_mfield() > 0) { InputStreamReader isr; String dataString = directory_name + "/calN.dat"; if (!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); isr = new InputStreamReader(response.getEntity().getContent()); } else { // Read from assets isr = new InputStreamReader(context.getAssets().open(dataString)); } BufferedReader reader = new BufferedReader(isr, buffer_size); String line = reader.readLine(); set_calN(Integer.parseInt(line)); reader.close(); } Log.d(DEBUG_TAG, "Finished reading calN.dat"); // Reading uL data if (get_Q_uL() > 0) { uL_vector = new Complex[get_Q_uL()][get_calN()]; for (int q_uL = 0; q_uL < get_Q_uL(); q_uL++) { BinaryReader dis; String dataString = directory_name + "/uL_" + String.format("%03d", q_uL) + ".bin"; if (!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); dis = new BinaryReader(response.getEntity().getContent()); } else { // Read from assets dis = new BinaryReader(context.getAssets().open(dataString)); } float[] Rdata = new float[get_calN()]; float[] Idata = new float[get_calN()]; Rdata = dis.ReadFloat(get_calN()); // read in the real part Idata = dis.ReadFloat(get_calN()); // read in the imagine part for (int i = 0; i < get_calN(); i++) uL_vector[q_uL][i] = new Complex((double) Rdata[i], (double) Idata[i]); dis.close(); } } Log.d(DEBUG_TAG, "Finished reading uL.dat"); // Read in Z data if (get_mfield() > 0) { Z_vector = new Complex[get_mfield()][get_n_basis_functions()][get_calN()]; for (int imf = 0; imf < get_mfield(); imf++) for (int inbfs = 0; inbfs < get_n_basis_functions(); inbfs++) { BinaryReader dis; String dataString = directory_name + "/Z_" + String.format("%03d", imf) + "_" + String.format("%03d", inbfs) + ".bin"; if (!isAssetFile) { HttpGet request = new HttpGet(dataString); HttpResponse response = client.execute(request); dis = new BinaryReader(response.getEntity().getContent()); } else { // Read from assets dis = new BinaryReader(context.getAssets().open(dataString)); } float[] Rdata = new float[get_calN()]; float[] Idata = new float[get_calN()]; Rdata = dis.ReadFloat(get_calN()); // read in the real part Idata = dis.ReadFloat(get_calN()); // read in the imagine part for (int i = 0; i < get_calN(); i++) Z_vector[imf][inbfs][i] = new Complex((double) Rdata[i], (double) Idata[i]); dis.close(); } } Log.d(DEBUG_TAG, "Finished reading Z.dat"); initialize_data_vectors(); }
From source file:rb.app.ComplexRBSystem.java
@Override public double RB_solve(int N) { current_N = N;//from w w w .j a v a2s . com 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 void loadOfflineData_rbappmit(AModelManager m) throws IOException { isReal = false;// ww w. j a v a2s. c o m // Find out dimension of the RB space { BufferedReader reader = m.getBufReader("n_bfs.dat"); set_n_basis_functions(Integer.parseInt(reader.readLine())); reader.close(); } Log.d(DEBUG_TAG, "Finished reading n_bfs.dat"); // Read in output data if (getNumOutputs() > 0) { // Get output dual norms { RB_output_vectors = (FieldVector<Complex>[][]) new FieldVector<?>[getNumOutputs()][]; output_dual_norms = new Complex[getNumOutputs()][]; for (int i = 0; i < getNumOutputs(); i++) { String[] dual_norms_tokens; { BufferedReader reader = m .getBufReader("output_" + String.format("%03d", i) + "_dual_norms.dat"); dual_norms_tokens = reader.readLine().split(" "); reader.close(); } { int Q_l_hat = getQl(i) * (getQl(i) + 1) / 2; output_dual_norms[i] = new Complex[Q_l_hat]; for (int q = 0; q < Q_l_hat; q++) { output_dual_norms[i][q] = new Complex(Double.parseDouble(dual_norms_tokens[q]), Double.parseDouble(dual_norms_tokens[Q_l_hat + q])); } } { RB_output_vectors[i] = (FieldVector<Complex>[]) new FieldVector<?>[getQl(i)]; String[] output_i_tokens; for (int q_l = 0; q_l < getQl(i); q_l++) { // Now read in the RB output vectors { BufferedReader reader_i = m.getBufReader("output_" + String.format("%03d", i) + "_" + String.format("%03d", q_l) + ".dat"); output_i_tokens = reader_i.readLine().split(" "); reader_i.close(); } RB_output_vectors[i][q_l] = new ArrayFieldVector<Complex>(getNBF(), new Complex(0d, 0d)); for (int j = 0; j < getNBF(); j++) { RB_output_vectors[i][q_l].setEntry(j, new Complex(Double.parseDouble(output_i_tokens[j]), Double.parseDouble(output_i_tokens[getNBF() + j]))); } } } } } Log.d(DEBUG_TAG, "Finished reading output data"); } else Log.d(DEBUG_TAG, "No output data set. (get_n_outputs() == 0)"); /* * // Read in the inner product matrix { InputStreamReader isr; String * dataString = directory_name + "/RB_inner_product_matrix.dat"; * * if(!isAssetFile) { HttpGet request = new HttpGet(dataString); * HttpResponse response = client.execute(request); isr = new * InputStreamReader(response.getEntity() .getContent()); } else { // * Read from assets isr = new InputStreamReader( * context.getAssets().open(dataString)); } BufferedReader * BufferedReader reader = new BufferedReader(isr,buffer_size); * * String String line = reader.readLine(); String[] tokens = * line.split(" "); * * // Set the size of the inner product matrix RB_inner_product_matrix = * new Array2DRowRealMatrix(n_bfs, n_bfs); * * // Fill the matrix int count = 0; for (int i = 0; i < n_bfs; i++) for * (int j = 0; j < n_bfs; j++) { RB_inner_product_matrix.setEntry(i, j, * Double .parseDouble(tokens[count])); count++; } reader.close(); } * * Log.d(DEBUG_TAG, "Finished reading RB_inner_product_matrix.dat"); */ // Read in the F_q vectors { RB_F_q_vector = (FieldVector<Complex>[]) new FieldVector<?>[getQf()]; String[] tokens; for (int q_f = 0; q_f < getQf(); q_f++) { { BufferedReader reader = m.getBufReader("RB_F_" + String.format("%03d", q_f) + ".dat"); tokens = reader.readLine().split(" "); reader.close(); } // Set the size of the inner product matrix RB_F_q_vector[q_f] = new ArrayFieldVector<Complex>(getNBF(), new Complex(0d, 0d)); // Fill the vector for (int i = 0; i < getNBF(); i++) { RB_F_q_vector[q_f].setEntry(i, new Complex(Double.parseDouble(tokens[i]), Double.parseDouble(tokens[getNBF() + i]))); } } Log.d(DEBUG_TAG, "Finished reading RB_F_q data"); } // Read in the A_q matrices { RB_A_q_vector = (FieldMatrix<Complex>[]) Array.newInstance(FieldMatrix.class, getQa());// (FieldMatrix<Complex>[]) // new // FieldMatrix<?>[get_Q_a()]; String[] tokens; for (int q_a = 0; q_a < getQa(); q_a++) { { BufferedReader reader = m.getBufReader("RB_A_" + String.format("%03d", q_a) + ".dat"); tokens = reader.readLine().split(" "); reader.close(); } // Set the size of the inner product matrix RB_A_q_vector[q_a] = new Array2DRowFieldMatrix<Complex>((new Complex(0, 0)).getField(), getNBF(), getNBF()); // Fill the vector int count = 0; for (int i = 0; i < getNBF(); i++) for (int j = 0; j < getNBF(); j++) { RB_A_q_vector[q_a].setEntry(i, j, new Complex(Double.parseDouble(tokens[count]), Double.parseDouble(tokens[count + getNBF() * getNBF()]))); count++; } } Log.d(DEBUG_TAG, "Finished reading RB_A_q data"); } // Read in F_q representor norm data { BufferedReader reader = m.getBufReader("Fq_norms.dat"); String[] tokens = reader.readLine().split(" "); reader.close(); // Declare the array int Q_f_hat = getQf() * (getQf() + 1) / 2; Fq_representor_norms = new Complex[Q_f_hat]; // Fill it for (int i = 0; i < Q_f_hat; i++) { Fq_representor_norms[i] = new Complex(Double.parseDouble(tokens[i * 2 + 0]), Double.parseDouble(tokens[i * 2 + 1])); } Log.d(DEBUG_TAG, "Finished reading Fq_norms.dat"); } // Read in Fq_Aq representor norm data { BufferedReader reader = m.getBufReader("Fq_Aq_norms.dat"); String[] tokens = reader.readLine().split(" "); reader.close(); reader = null; // Declare the array Fq_Aq_representor_norms = new Complex[getQf()][getQa()][getNBF()]; double[][][] Rdata = new double[getQf()][getQa()][getNBF()]; double[][][] Idata = new double[getQf()][getQa()][getNBF()]; // Fill it int count = 0; for (int q_f = 0; q_f < getQf(); q_f++) for (int q_a = 0; q_a < getQa(); q_a++) { for (int i = 0; i < getNBF(); i++) { Rdata[q_f][q_a][i] = Double.parseDouble(tokens[count]); count++; } for (int i = 0; i < getNBF(); i++) { Idata[q_f][q_a][i] = Double.parseDouble(tokens[count]); count++; } } for (int q_f = 0; q_f < getQf(); q_f++) for (int q_a = 0; q_a < getQa(); q_a++) for (int i = 0; i < getNBF(); i++) Fq_Aq_representor_norms[q_f][q_a][i] = new Complex(Rdata[q_f][q_a][i], Idata[q_f][q_a][i]); Log.d(DEBUG_TAG, "Finished reading Fq_Aq_norms.dat"); } MathObjectReader mr = m.getMathObjReader(); // Read in Aq_Aq representor norm data { // Declare the array int Q_a_hat = getQa() * (getQa() + 1) / 2; Aq_Aq_representor_norms = new Complex[Q_a_hat][getNBF()][getNBF()]; int count = 0; double[][] Rdata2 = null, Idata2 = null; for (int i = 0; i < getQa(); i++) for (int j = i; j < getQa(); j++) { String file = "Aq_Aq_" + String.format("%03d", i) + "_" + String.format("%03d", j) + "_norms.bin"; int n = getNBF(); InputStream in = m.getInStream(file); try { Rdata2 = mr.readRawDoubleMatrix(in, n, n); Idata2 = mr.readRawDoubleMatrix(in, n, n); } catch (IOException io) { Log.e("ComplexRBSystem", "IOException with file " + file, io); } finally { in.close(); in = null; } for (int k = 0; k < n; k++) for (int l = 0; l < n; l++) Aq_Aq_representor_norms[count][k][l] = new Complex(Rdata2[k][l], Idata2[k][l]); count++; } Log.d(DEBUG_TAG, "Finished reading Aq_Aq_norms.dat"); } // // Read calN number // { // if (getNumFields() > 0) { // BufferedReader reader = m.getBufReader("calN.dat"); // // String line = reader.readLine(); // // set_calN(Integer.parseInt(line)); // reader.close(); // } // // Log.d(DEBUG_TAG, "Finished reading calN.dat"); // } int n = getGeometry().getNumVertices(); // Reading uL data { if (getQuL() > 0) { uL_vector = new Complex[getQuL()][n]; float[] Rdata3, Idata3; for (int q_uL = 0; q_uL < getQuL(); q_uL++) { InputStream in = m.getInStream("uL_" + String.format("%03d", q_uL) + ".bin"); try { Rdata3 = mr.readRawFloatVector(in, n); Idata3 = mr.readRawFloatVector(in, n); } finally { in.close(); } for (int i = 0; i < n; i++) uL_vector[q_uL][i] = new Complex(Rdata3[i], Idata3[i]); } } Log.d(DEBUG_TAG, "Finished reading uL.dat"); } // Read in Z data { if (getNumDoFFields() > 0) { Z_vector = new Complex[getNumDoFFields()][getNBF()][n]; float[] Rdata3, Idata3; for (int imf = 0; imf < getNumDoFFields(); imf++) for (int inbfs = 0; inbfs < getNBF(); inbfs++) { InputStream in = m.getInStream( "Z_" + String.format("%03d", imf) + "_" + String.format("%03d", inbfs) + ".bin"); try { Rdata3 = mr.readRawFloatVector(in, n); Idata3 = mr.readRawFloatVector(in, n); } finally { in.close(); } for (int i = 0; i < n; i++) Z_vector[imf][inbfs][i] = new Complex(Rdata3[i], Idata3[i]); } } Log.d(DEBUG_TAG, "Finished reading Z.dat"); } initialize_data_vectors(); }
From source file:rb.ComplexRBSystem.java
@Override public double RB_solve(int N) { current_N = N;//w w w. j a v a 2 s . 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); }