Commit 184bdf09 authored by Daniel Brown's avatar Daniel Brown
Browse files

Fixing mistake in detecting whether a nodes q value has changed or not. Was...

Fixing mistake in detecting whether a nodes q value has changed or not. Was leading to a recomputation of k matrices when not needed. Reduced an aLIGO file runtime from 23s to 12s though. Looked further at sparse vector input to klu, turns out I read the code wrong before and not as straight forward so removing the usage from quantum calculation, shouldn't affect results though.
parent 0fd71548
This diff is collapsed.
......@@ -548,16 +548,16 @@ Int KLU_tsolve_sparse
for (nb = 0; nb < nBblocks; nb++)
{
Int start = Bblocks[2*nb];
Int end = Bblocks[2*nb+1];
Int num = Bblocks[2*nb+1];
for (chunk = start ; chunk < end ; chunk += 4)
for (chunk = start ; chunk < start + num ; chunk += 4)
{
/* ------------------------------------------------------------------ */
/* get the size of the current chunk */
/* ------------------------------------------------------------------ */
nr = MIN (end - chunk, 4) ;
nr = MIN (start + num - chunk, 4) ;
/* ------------------------------------------------------------------ */
/* permute the right hand side, X = Q'*B */
......
......@@ -994,6 +994,7 @@ int main(int argc, char *argv[]) {
build_ccs_matrix(&M_ifo_car);
get_ccs_matrix(&M_ifo_car);
set_light_input_ports();
//get_ccs_matrix_sparse_input_blocks();
if (M_ifo_sig.num_frequencies > 0) {
// The signal matrix also depends on the number of signal frequencies and also the number
......
......@@ -384,6 +384,8 @@ typedef struct node {
bool q_changed; //!< True if the q value has been changed to what it previously was
complex_t qx; //!< x component of beam q parameter
complex_t qy; //!< y component of beam q parameter
complex_t _prev_qx; //!< x component of beam q parameter before retrace
complex_t _prev_qy; //!< y component of beam q parameter before retrace
int component_index; //!< component index
double *n; //!< refractive index of an attached space
bool q_is_set; //!< is Gaussian beam parameter set?
......
......@@ -3006,7 +3006,9 @@ void set_k_all_components(bool rebuilding) {
int mm = 0;
int mbs = 0;
int ms = 0;
int ml = 0;
int ml = 0;
int mmd = 0;
int md = 0;
initonly = noinit = 0;
......@@ -3028,9 +3030,9 @@ void set_k_all_components(bool rebuilding) {
for (i = 0; i < inter.num_mirrors; i++) {
mirror_t *mirror = &(inter.mirror_list[i]);
if (!rebuilding || mirror->rebuild)
if (!rebuilding || mirror->rebuild) {
mm = set_k_mirror(i);
else{
} else {
// no known rebuilding, still need to check if the q values have
// changed
node_t n1 = inter.node_list[mirror->node1_index];
......@@ -3040,17 +3042,18 @@ void set_k_all_components(bool rebuilding) {
//warn("mirror %s node %s changed %d\n",mirror->name, n2.name, n2.q_changed);
if(((!n1.gnd_node) && (n1.q_changed)) ||
((!n2.gnd_node) && (n2.q_changed)))
((!n2.gnd_node) && (n2.q_changed))) {
mm = set_k_mirror(i);
}
}
}
for (i = 0; i < inter.num_beamsplitters; i++) {
beamsplitter_t *bs = &(inter.bs_list[i]);
if (!rebuilding || bs->rebuild)
if (!rebuilding || bs->rebuild){
mbs = set_k_beamsplitter(i);
else{
}else{
// no known rebuilding, still need to check if the q values have
// changed
node_t n1 = inter.node_list[bs->node1_index];
......@@ -3061,24 +3064,26 @@ void set_k_all_components(bool rebuilding) {
if(((!n1.gnd_node) && (n1.q_changed))
|| ((!n2.gnd_node) && (n2.q_changed))
|| ((!n3.gnd_node) && (n3.q_changed))
|| ((!n4.gnd_node) && (n4.q_changed)))
mm = set_k_beamsplitter(i);
|| ((!n4.gnd_node) && (n4.q_changed))){
mbs = set_k_beamsplitter(i);
}
}
}
for (i = 0; i < inter.num_spaces; i++) {
space_t *s = &inter.space_list[i];
if (!rebuilding || s->rebuild)
if (!rebuilding || s->rebuild) {
ms = set_k_space(i);
else{
}else{
// no known rebuilding, still need to check if the q values have
// changed
node_t n1 = inter.node_list[s->node1_index];
node_t n2 = inter.node_list[s->node2_index];
if(((!n1.gnd_node) && (n1.q_changed))
|| ((!n2.gnd_node) && (n2.q_changed)))
mm = set_k_space(i);
|| ((!n2.gnd_node) && (n2.q_changed))){
ms = set_k_space(i);
}
}
if (inter.trace & 4 && ms) {
......@@ -3089,51 +3094,54 @@ void set_k_all_components(bool rebuilding) {
for (i = 0; i < inter.num_lenses; i++) {
lens_t *l = &inter.lens_list[i];
if (!rebuilding || l->rebuild)
if (!rebuilding || l->rebuild){
ml = set_k_lens(i);
else{
}else{
// no known rebuilding, still need to check if the q values have
// changed
node_t n1 = inter.node_list[l->node1_index];
node_t n2 = inter.node_list[l->node2_index];
if(((!n1.gnd_node) && (n1.q_changed))
|| ((!n2.gnd_node) && (n2.q_changed)))
mm = set_k_lens(i);
|| ((!n2.gnd_node) && (n2.q_changed))){
ml = set_k_lens(i);
}
}
}
for (i = 0; i < inter.num_modulators; i++) {
modulator_t *mod = &inter.modulator_list[i];
if (!rebuilding || mod->rebuild)
ml = set_k_modulator(i);
else{
if (!rebuilding || mod->rebuild){
mmd = set_k_modulator(i);
}else{
// no known rebuilding, still need to check if the q values have
// changed
node_t n1 = inter.node_list[mod->node1_index];
node_t n2 = inter.node_list[mod->node2_index];
if(((!n1.gnd_node) && (n1.q_changed))
|| ((!n2.gnd_node) && (n2.q_changed)))
mm = set_k_modulator(i);
|| ((!n2.gnd_node) && (n2.q_changed))){
mmd = set_k_modulator(i);
}
}
}
for (i = 0; i < inter.num_diodes; i++) {
diode_t *d = &inter.diode_list[i];
if (!rebuilding || d->rebuild)
ml = set_k_diode(i);
else{
if (!rebuilding || d->rebuild){
md = set_k_diode(i);
}else{
// no known rebuilding, still need to check if the q values have
// changed
node_t n1 = inter.node_list[d->node1_index];
node_t n2 = inter.node_list[d->node2_index];
if(((!n1.gnd_node) && (n1.q_changed))
|| ((!n2.gnd_node) && (n2.q_changed)))
mm = set_k_diode(i);
|| ((!n2.gnd_node) && (n2.q_changed))){
md = set_k_diode(i);
}
}
}
......@@ -3155,7 +3163,8 @@ void set_k_all_components(bool rebuilding) {
(void) mbs;
(void) ms;
(void) ml;
(void) mmd;
(void) md;
}
//! Set up an ABCD mirror
......@@ -5428,7 +5437,6 @@ void set_q(int node_index, int component_index, complex_t qx, complex_t qy) {
message("--- setting q for node %s\n", get_node_name(node_index));
}
node->q_changed = !ceq(qx, node->qx) || !ceq(qy, node->qy);
node->qx = qx;
node->qy = qy;
node->component_index = component_index;
......@@ -5485,7 +5493,6 @@ void set_q(int node_index, int component_index, complex_t qx, complex_t qy) {
message(" overwriting second q value at space %s\n", space.name);
}
node->q_changed = !ceq(qxt, node->qx) || !ceq(qyt, node->qy);
node->qx = qxt;
node->qy = qyt;
node->component_index = component1_index;
......@@ -5501,7 +5508,6 @@ void set_q(int node_index, int component_index, complex_t qx, complex_t qy) {
message(" setting second q value at space %s\n", space.name);
}
node->q_changed = !ceq(qxt, node->qx) || !ceq(qyt, node->qy);
node->qx = qxt;
node->qy = qyt;
node->component_index = component1_index;
......@@ -5564,7 +5570,6 @@ void set_q(int node_index, int component_index, complex_t qx, complex_t qy) {
message(" overwriting second q value at space %s\n", space.name);
}
node->q_changed = !ceq(qxt, node->qx) || !ceq(qyt, node->qy);
node->qx = qxt;
node->qy = qyt;
node->component_index = component2_index;
......@@ -5580,7 +5585,6 @@ void set_q(int node_index, int component_index, complex_t qx, complex_t qy) {
message(" setting second q value at space %s\n", space.name);
}
node->q_changed = !ceq(qxt, node->qx) || !ceq(qyt, node->qy);
node->qx = qxt;
node->qy = qyt;
node->component_index = component2_index;
......
......@@ -1583,6 +1583,7 @@ void data_point_new(){
}
// re fill the matrix with the latest values for the data point
fill_ccs_matrix(&M_ifo_car);
factor_matrix(STANDARD);
clear_rhs(STANDARD);
......@@ -1656,7 +1657,11 @@ void data_point_new(){
complex_t input_amplitude = z_by_ph(co(sqrt(2*light_input->I0 / init.epsilon_0_c), 0.0), light_input->phase);
input_amplitude = z_by_z(input_amplitude, light_input->power_coeff_list[j]);
set_rhs(rhs_idx_car + get_node_rhs_idx(1, light_input->f->index, j, M_ifo_car.num_frequencies), input_amplitude, STANDARD);
if(input_amplitude.im == 0.0 && input_amplitude.re == 0.0)
continue;
int idx = rhs_idx_car + get_node_rhs_idx(1, light_input->f->index, j, M_ifo_car.num_frequencies);
set_rhs(idx, input_amplitude, STANDARD);
// if there is a signal on this laser then add in the sideband source in RHS
if(signal){
......@@ -3478,7 +3483,24 @@ void rebuild_all(void) {
// retrace all
if (inter.retrace) {
_tid = startTimer("REBUILD_RETRACE");
for (i = 0; i < inter.num_nodes; i++) {
inter.node_list[i]._prev_qx = inter.node_list[i].qx;
inter.node_list[i]._prev_qy = inter.node_list[i].qy;
}
retrace();
for (i = 0; i < inter.num_nodes; i++) {
if(!ceq(inter.node_list[i]._prev_qx, inter.node_list[i].qx) ||
!ceq(inter.node_list[i]._prev_qy, inter.node_list[i].qy)) {
inter.node_list[i].q_changed = true;
} else {
inter.node_list[i].q_changed = false;
}
}
endTimer(_tid);
}
}
......
......@@ -4294,6 +4294,71 @@ void build_ccs_matrix(ifo_matrix_vars_t *matrix){
}
}
/**
* Typically RHS vector is also sparse for Finesse models, thus we can be efficient
* and determine which values are non-zero beforehand.
*
* @param matrix
*/
void get_ccs_matrix_sparse_input_blocks() {
M_ifo_car.input_blocks = (long*)malloc(sizeof(long) * inter.num_light_inputs * inter.num_fields);
M_ifo_car.input_num_blocks = 0;
M_ifo_sig.input_num_blocks = 0;
int i, j;
for(i=0; i<inter.num_light_inputs; i++){
light_in_t *light_input = &(inter.light_in_list[i]);
node_t* n = &inter.node_list[light_input->node_index];
int rhs_idx_car = 0;
// first get the
if(light_input->node_port == 2) { // first port in the input
rhs_idx_car = M_ifo_car.node_rhs_idx_2[n->list_index];
} else if(light_input->node_port == 1) { // second port is the input
rhs_idx_car = M_ifo_car.node_rhs_idx_1[n->list_index];
} else
bug_error("Node did not have direction set when filling in RHS vector");
int inblock = false, idx = 0;
for (j = 0; j < inter.num_fields; j++) {
idx = rhs_idx_car + get_port_rhs_idx(light_input->f->index, j);
complex_t v = light_input->power_coeff_list[j];
if(v.re == 0.0 && v.im == 0.0 && inblock) {
// we're in a block but now there is no power
inblock = false;
} else if((v.re != 0.0 || v.im != 0.0) && !inblock) {
// have a non-zero value and not in block so create one
M_ifo_car.input_num_blocks++;
M_ifo_car.input_blocks[2*(M_ifo_car.input_num_blocks-1)] = idx;
M_ifo_car.input_blocks[2*(M_ifo_car.input_num_blocks-1) + 1] = 1;
inblock = true;
} else if(inblock) {
M_ifo_car.input_blocks[2*(M_ifo_car.input_num_blocks-1) + 1] += 1;
}
}
}
message("\n");
warn("-----------------------------\n");
warn("- Carrier RHS vector blocks\n");
warn("-----------------------------\n");
for(i=0; i<M_ifo_car.input_num_blocks; i++) {
warn("Block %i: %i -> %i\n", i, M_ifo_car.input_blocks[2*i], M_ifo_car.input_blocks[2*i+1]);
}
warn("-----------------------------\n");
}
void get_ccs_matrix(ifo_matrix_vars_t *matrix){
// reset pseudo-column linked list
pseudocolumn_t *curr = matrix->pcol_head;
......@@ -4589,14 +4654,18 @@ void solve_ccs_matrix(ifo_matrix_vars_t *matrix, double *rhs, bool transpose, bo
if(options.sparse_solver == KLU_FULL){
klu_objs_t *ko = (klu_objs_t*)matrix->solver_opts;
if(transpose)
if(conjugate)
klu_zl_tsolve(ko->Symbolic, ko->Numeric, matrix->M.num_eqns, 1, rhs, 1, &ko->Common);
else
klu_zl_tsolve(ko->Symbolic, ko->Numeric, matrix->M.num_eqns, 1, rhs, 0, &ko->Common);
else
if(transpose) {
//if (!options.use_coupling_reduction)
klu_zl_tsolve(ko->Symbolic, ko->Numeric, matrix->M.num_eqns, 1, rhs, conjugate, &ko->Common);
//else
// klu_zl_tsolve_sparse(ko->Symbolic, ko->Numeric, matrix->M.num_eqns, 1, rhs, matrix->input_num_blocks, matrix->input_blocks, conjugate, &ko->Common);
} else {
//if (!options.use_coupling_reduction)
klu_zl_solve(ko->Symbolic, ko->Numeric, matrix->M.num_eqns, 1, rhs, &ko->Common);
//else
//klu_zl_solve_sparse(ko->Symbolic, ko->Numeric, matrix->M.num_eqns, 1, rhs, matrix->input_num_blocks, matrix->input_blocks, &ko->Common);
}
if(ko->Common.status != 0)
gerror("KLU error whilst solving: %i", ko->Common.status);
#if INCLUDE_NICSLU == 1
......
......@@ -89,9 +89,13 @@ typedef struct ifo_matrix_vars{
uint64_t *output_rhs_idx; /** For a given output index the RHS index of the port this output is attached to is stored. This must be 64 bit as for homodyne detectors we store 2 ints in one 64 bit type. */
int *block_rhs_idx;
int *feedback_rhs_idx;
int input_num_blocks; /** Number of blocks of non-zero elements in the input array */
long *input_blocks; /** Size 2*input_num_blocks, n=start index of non-zeros, n+1=end index of non-zeros */
} ifo_matrix_vars_t;
void get_ccs_matrix_sparse_input_blocks();
void dump_matrix_ccs(matrix_ccs_t *matrix, const char* fname);
void dump_matrix_ccs_labels(ifo_matrix_vars_t *ifo_matrix, const char* fname);
......
......@@ -1446,25 +1446,26 @@ void compute_squeezing_factor(light_out_t *sd, double *r, double *phi, double *d
memset(quant_w, 0, sizeof(complex_t)*M_ifo_sig.M.num_eqns);
memset(quant_v, 0, sizeof(complex_t)*M_ifo_sig.M.num_eqns);
long nblocks = 1;
long blocks[2] = {0};
// fill selection vector to select column
quant_w[idx[j]] = complex_1;
// store which elements of the input matrix are non-zero
blocks[0] = idx[j];
blocks[1] = idx[j]+1;
// long nblocks = 1;
// long blocks[2] = {0};
// // store which elements of the input vector are non-zero
// blocks[0] = idx[j]; // start index
// blocks[1] = 2; // number of elements
// compute the w = M^-1^T s, to get the weighting vector
solve_ccs_matrix_sparse_input(&M_ifo_sig, (double*)quant_w, true, nblocks, blocks);
solve_ccs_matrix(&M_ifo_sig, (double*)quant_w, true, false);
// apply input noise matrix the weighting vectors
// i = M_qn * w
ccs_zgemv(&M_qni, quant_w, quant_v);
// now compute the final vector <b b^T>s = M^-1 v
solve_ccs_matrix_sparse_input(&M_ifo_sig, (double*)quant_v, false, nblocks, blocks);
solve_ccs_matrix(&M_ifo_sig, (double*)quant_v, false, false);
b[j][0] = z_by_x(quant_v[idx[0]],2);
b[j][1] = z_by_x(quant_v[idx[1]],2);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment