text:
#include <gsl/gsl_vector.h>
#include <stdio.h>
#include "vector.h"
#include "api.h"

 static void ke_vector_alloc(sml_t* sml) {
	int n = sml_pop_int(sml);
	gsl_vector * v = gsl_vector_alloc(n);
	sml_mem_inc(sml);
	sml_push_vector(sml, v)
}

 static void ke_vector_int_alloc(sml_t* sml) {
	int n = sml_pop_int(sml);
	gsl_vector_int * v = gsl_vector_int_alloc(n);
	sml_mem_inc(sml);
	sml_push_vector_int(sml, v);
}

 static void ke_vector_free(sml_t* sml) {
	token_t * tokp = sml_pop_token(sml);
	void * v = sml_get_ptr(tokp);
	sml_free_ptr(sml, v);
	sml_set_ptr_null(sml,tokp);
}

 static void ke_vector_int_free(sml_t* sml) {
	 token_t * tokp = sml_pop_token(sml);
	 void * v = sml_get_ptr(tokp);
	 sml_free_ptr(sml, v);
	 sml_set_ptr_null(sml,tokp);
}

void ke_vector_get(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	int i = sml_pop_int(sml);
	double r = gsl_vector_get(v, i);
	sml_push_real(sml, r);
}

void ke_vector_int_get(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	int idx = sml_pop_int(sml);
	int i = gsl_vector_int_get(v, idx);
	sml_push_int(sml, i);
}

void ke_vector_prop_set(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_vector * v = sml_pop_vector(sml);
	int i = sml_pop_int(sml);
	gsl_vector_set(v, i, x);
}

void ke_vector_int_set(sml_t* sml) {
	int x = sml_pop_int(sml);
	int i = sml_pop_int(sml);
	gsl_vector_int * v = sml_pop_vector_int(sml);
	gsl_vector_int_set(v, i, x);
}

void ke_vector_prop_int_set(sml_t* sml) {
	int x = sml_pop_int(sml);
	gsl_vector_int * v = sml_pop_vector_int(sml);
	int i = sml_pop_int(sml);
	gsl_vector_int_set(v,i,x);
}

void ke_vector_set(sml_t* sml) {
	double x = sml_pop_real(sml);
	int i = sml_pop_int(sml);
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector_set(v,i,x);
}

void ke_vector_put(sml_t* sml) {
	int top = sml_get_top(sml);
	int n = sml_get_args(sml);
	gsl_vector * v = sml_peek_vector(sml, (top - n));
	for (size_t i = 0; i < n - 1; ++i) {
		if (i == v->size) break;
		double x = sml_peek_real(sml, (top - n + i + 1));
		gsl_vector_set(v, i, x);
	}
	top = top - n;
	sml_set_top(sml, top);
}

static void ke_vector_int_put(sml_t* sml) {
	int top = sml_get_top(sml);
	int n = sml_get_args(sml);
	gsl_vector_int * v = sml_peek_vector_int(sml, (top - n));
	for (size_t i = 0; i < n - 1; ++i) {
		if (i == v->size) break;
		int x = sml_peek_int(sml, (top - n + i + 1));
		gsl_vector_int_set(v, i, x);
	}
	top = top - n;
	sml_set_top(sml, top);
}

 static void ke_vector_set_all(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector_set_all(v,x);
}

 static void ke_vector_int_set_all(sml_t* sml) {
	 int x = sml_pop_int(sml);
	 gsl_vector_int * v = sml_pop_vector_int(sml);
	 gsl_vector_int_set_all(v, x);
}

 static void ke_vector_set_zero(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector_set_zero(v);
}

 static void ke_vector_int_set_zero(sml_t* sml) {
	 gsl_vector_int * v = sml_pop_vector_int(sml);
	 gsl_vector_int_set_zero(v);
}

 static void ke_vector_set_basis(sml_t* sml) {
	int i = sml_pop_int(sml);
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector_set_basis(v, i);
}

 static void ke_vector_int_set_basis(sml_t* sml) {
	 int i = sml_pop_int(sml);
	 gsl_vector_int * v = sml_pop_vector_int(sml);
	 gsl_vector_int_set_basis(v, i);
}

 static void ke_vector_add(sml_t* sml) {
	gsl_vector * b = sml_pop_vector(sml);
	gsl_vector * a = sml_pop_vector(sml);
	gsl_vector_add(a,b);
}

 static void ke_vector_int_add(sml_t* sml) {
	gsl_vector_int * b = sml_pop_vector_int(sml);
	gsl_vector_int * a = sml_pop_vector_int(sml);
	gsl_vector_int_add(a, b);
}

 static void ke_vector_sub(sml_t* sml) {
	gsl_vector * b = sml_pop_vector(sml);
	gsl_vector * a = sml_pop_vector(sml);
	gsl_vector_sub(a, b);
}

 static void ke_vector_int_sub(sml_t* sml) {
	gsl_vector_int * b = sml_pop_vector_int(sml);
	gsl_vector_int * a = sml_pop_vector_int(sml);
	gsl_vector_int_sub(a, b);
}

 static void ke_vector_mul(sml_t* sml) {
	gsl_vector * b = sml_pop_vector(sml);
	gsl_vector * a = sml_pop_vector(sml);
	gsl_vector_mul(a, b);
}

 static void ke_vector_int_mul(sml_t* sml) {
	gsl_vector_int * b = sml_pop_vector_int(sml);
	gsl_vector_int * a = sml_pop_vector_int(sml);
	gsl_vector_int_mul(a, b);
 }

 static void ke_vector_div(sml_t* sml) {
	gsl_vector * b = sml_pop_vector(sml);
	gsl_vector * a = sml_pop_vector(sml);
	gsl_vector_div(a, b);
}

 static void ke_vector_int_div(sml_t* sml) {
	gsl_vector_int * b = sml_pop_vector_int(sml);
	gsl_vector_int * a = sml_pop_vector_int(sml);
	gsl_vector_int_div(a, b);
}

 static void ke_vector_scale(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector_scale(v,x);
}

 static void ke_vector_int_scale(sml_t* sml) {
	 double x = sml_pop_real(sml);
	 gsl_vector_int * v = sml_pop_vector_int(sml);
	 gsl_vector_int_scale(v, x);
}

 static void ke_vector_add_constant(sml_t* sml) {
	double x = sml_pop_real(sml);
	gsl_vector * a = sml_pop_vector(sml);
	gsl_vector_add_constant(a,x);
}

 static void ke_vector_int_add_constant(sml_t* sml) {
	 double x = sml_pop_real(sml);
	 gsl_vector_int * a = sml_pop_vector_int(sml);
	 gsl_vector_int_add_constant(a, x);
}

 static void ke_vector_reverse(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector_reverse(v);
}

 static void ke_vector_int_reverse(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	gsl_vector_int_reverse(v);
}

 static void ke_vector_swap_elements(sml_t* sml) {
	int j = sml_pop_int(sml);
	int i = sml_pop_int(sml);
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector_swap_elements(v,i,j);
}

 static void ke_vector_int_swap_elements(sml_t* sml) {
	 int j = sml_pop_int(sml);
	 int i = sml_pop_int(sml);
	 gsl_vector_int * v = sml_pop_vector_int(sml);
	 gsl_vector_int_swap_elements(v, i, j);
}

 static void ke_vector_memcpy(sml_t* sml) {
	gsl_vector * src = sml_pop_vector(sml);
	gsl_vector * dest = sml_pop_vector(sml);
	gsl_vector_memcpy(dest, src);
}

 static void ke_vector_int_memcpy(sml_t* sml) {
	gsl_vector_int * src = sml_pop_vector_int(sml);
	gsl_vector_int * dest = sml_pop_vector_int(sml);
	gsl_vector_int_memcpy(dest, src);
}

 static void ke_vector_swap(sml_t* sml) {
	gsl_vector * w = sml_pop_vector(sml);
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector_swap(v,w);
}

 static void ke_vector_int_swap(sml_t* sml) {
	gsl_vector_int * w = sml_pop_vector_int(sml);
	gsl_vector_int * v = sml_pop_vector_int(sml);
	gsl_vector_int_swap(v, w);
}

 static void ke_vector_min(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	double r = gsl_vector_min(v);
	sml_push_real(sml,r);
}

 static void ke_vector_int_min(sml_t* sml) {
	 gsl_vector_int * v = sml_pop_vector_int(sml);
	 int i = gsl_vector_int_min(v);
	 sml_push_int(sml, i);
}

 static void ke_vector_max(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	double r = gsl_vector_max(v);
	sml_push_real(sml, r);
}

 static void ke_vector_int_max(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	int i = gsl_vector_int_max(v);
	sml_push_int(sml, i);
}

 static void ke_vector_isnull(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	int i = gsl_vector_isnull(v);
	sml_push_int(sml, i);
}

 static void ke_vector_int_isnull(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	int i = gsl_vector_int_isnull(v);
	sml_push_int(sml, i);
}

 static void ke_vector_ispos(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	int i = gsl_vector_ispos(v);
	sml_push_int(sml, i);
}

 static void ke_vector_int_ispos(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	int i = gsl_vector_int_ispos(v);
	sml_push_int(sml, i);
}

 static void ke_vector_isneg(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	int i = gsl_vector_isneg(v);
	sml_push_int(sml, i);
}

 static void ke_vector_int_isneg(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	int i = gsl_vector_int_isneg(v);
	sml_push_int(sml, i);
}

 static void ke_vector_isnonneg(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	int i = gsl_vector_isnonneg(v);
	sml_push_int(sml, i);
}

 static void ke_vector_int_isnonneg(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	int i = gsl_vector_int_isnonneg(v);
	sml_push_int(sml, i);
}

 static void ke_vector_equal(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	gsl_vector * u = sml_pop_vector(sml);
	int i  = gsl_vector_equal(u,v);
	sml_push_int(sml, i);
 }

 static void ke_vector_int_equal(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	gsl_vector_int * u = sml_pop_vector_int(sml);
	int i = gsl_vector_int_equal(u, v);
	sml_push_int(sml, i);
}

 static void ke_vector_fscanf(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	char * filename = sml_pop_str(sml);
    FILE * f = fopen(filename, "r");
    gsl_vector_fscanf(f, v);
    fclose(f);
}

 static void ke_vector_int_fscanf(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	char * filename = sml_pop_str(sml);
	FILE * f = fopen(filename, "r");
	gsl_vector_int_fscanf(f, v);
	fclose(f);
}

 static void ke_vector_fprintf(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	char * filename = sml_pop_str(sml);
	FILE * f = fopen(filename, "w");
    gsl_vector_fprintf(f, v,"%5g");
    fclose(f);
}

 static void ke_vector_int_fprintf(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	char * filename = sml_pop_str(sml);
	FILE * f = fopen(filename, "w");
	gsl_vector_int_fprintf(f, v, "%5g");
	fclose(f);
}

 static void ke_vector_fread(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	char * filename = sml_pop_str(sml);
	FILE * f = fopen(filename, "r");
	gsl_vector_fread(f, v);
	fclose(f);
}

 static void ke_vector_int_fread(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	char * filename = sml_pop_str(sml);
	FILE * f = fopen(filename, "r");
	gsl_vector_int_fread(f, v);
	fclose(f);
}

 static void ke_vector_fwrite(sml_t* sml) {
	gsl_vector * v = sml_pop_vector(sml);
	char * filename = sml_pop_str(sml);
	FILE * f = fopen(filename, "w");
    gsl_vector_fwrite(f, v);
    fclose(f);
}

 static void ke_vector_int_fwrite(sml_t* sml) {
	gsl_vector_int * v = sml_pop_vector_int(sml);
	char * filename = sml_pop_str(sml);
	FILE * f = fopen(filename, "w");
	gsl_vector_int_fwrite(f, v);
	fclose(f);
}

void ke_vector_hash(sml_t* sml) {
    ke_hash_add(sml, (fncp)&ke_vector_alloc, VECTOR);
    ke_hash_add(sml, (fncp)&ke_vector_alloc, VECTOR_ALLOC);
    ke_hash_add(sml, (fncp)&ke_vector_get, VECTOR_GET);
    ke_hash_add(sml, (fncp)&ke_vector_set, VECTOR_SET);
	ke_hash_add(sml, (fncp)&ke_vector_get, VECTOR_SGET);
	ke_hash_add(sml, (fncp)&ke_vector_set, VECTOR_SSET);
	ke_hash_add(sml, (fncp)&ke_vector_put, VECTOR_PUT);
	ke_hash_add(sml, (fncp)&ke_vector_free, VECTOR_FREE);
    ke_hash_add(sml, (fncp)&ke_vector_set_all, VECTOR_SET_ALL);
    ke_hash_add(sml, (fncp)&ke_vector_set_zero, VECTOR_SET_ZERO);
    ke_hash_add(sml, (fncp)&ke_vector_set_basis, VECTOR_SET_BASIS);
    ke_hash_add(sml, (fncp)&ke_vector_add, VECTOR_ADD);
    ke_hash_add(sml, (fncp)&ke_vector_sub, VECTOR_SUB);
    ke_hash_add(sml, (fncp)&ke_vector_mul, VECTOR_MUL);
    ke_hash_add(sml, (fncp)&ke_vector_div, VECTOR_DIV);
    ke_hash_add(sml, (fncp)&ke_vector_scale, VECTOR_SCALE);
    ke_hash_add(sml, (fncp)&ke_vector_add_constant, VECTOR_ADD_CONSTANT);
    ke_hash_add(sml, (fncp)&ke_vector_reverse, VECTOR_REVERSE);
    ke_hash_add(sml, (fncp)&ke_vector_swap_elements, VECTOR_SWAP_ELEMENTS);
    ke_hash_add(sml, (fncp)&ke_vector_memcpy, VECTOR_MEMCPY);
    ke_hash_add(sml, (fncp)&ke_vector_swap, VECTOR_SWAP);
    ke_hash_add(sml, (fncp)&ke_vector_min, VECTOR_MIN);
    ke_hash_add(sml, (fncp)&ke_vector_max, VECTOR_MAX);
    ke_hash_add(sml, (fncp)&ke_vector_isnull, VECTOR_ISNULL);
    ke_hash_add(sml, (fncp)&ke_vector_ispos, VECTOR_ISPOS);
    ke_hash_add(sml, (fncp)&ke_vector_isneg, VECTOR_ISNEG);
    ke_hash_add(sml, (fncp)&ke_vector_isnonneg, VECTOR_ISNONNEG);
    ke_hash_add(sml, (fncp)&ke_vector_equal, VECTOR_EQUAL);
    ke_hash_add(sml, (fncp)&ke_vector_fscanf, VECTOR_FSCANF);
    ke_hash_add(sml, (fncp)&ke_vector_fscanf, VECTOR_SCAN);
    ke_hash_add(sml, (fncp)&ke_vector_fprintf, VECTOR_FPRINTF);
    ke_hash_add(sml, (fncp)&ke_vector_fprintf, VECTOR_PRINT);
    ke_hash_add(sml, (fncp)&ke_vector_fread, VECTOR_FREAD);
    ke_hash_add(sml, (fncp)&ke_vector_fread, VECTOR_READ);
    ke_hash_add(sml, (fncp)&ke_vector_fwrite, VECTOR_FWRITE);
    ke_hash_add(sml, (fncp)&ke_vector_fwrite, VECTOR_WRITE);

	ke_hash_add(sml, (fncp)&ke_vector_int_alloc, VECTOR_INT);
	ke_hash_add(sml, (fncp)&ke_vector_int_alloc, VECTOR_INT_ALLOC);
	ke_hash_add(sml, (fncp)&ke_vector_int_get, VECTOR_INT_GET);
	ke_hash_add(sml, (fncp)&ke_vector_int_set, VECTOR_INT_SET);
	ke_hash_add(sml, (fncp)&ke_vector_int_put, VECTOR_INT_PUT);
	ke_hash_add(sml, (fncp)&ke_vector_int_free, VECTOR_INT_FREE);
	ke_hash_add(sml, (fncp)&ke_vector_int_set_all, VECTOR_INT_SET_ALL);
	ke_hash_add(sml, (fncp)&ke_vector_int_set_zero, VECTOR_INT_SET_ZERO);
	ke_hash_add(sml, (fncp)&ke_vector_int_set_basis, VECTOR_INT_SET_BASIS);
	ke_hash_add(sml, (fncp)&ke_vector_int_add, VECTOR_INT_ADD);
	ke_hash_add(sml, (fncp)&ke_vector_int_sub, VECTOR_INT_SUB);
	ke_hash_add(sml, (fncp)&ke_vector_int_mul, VECTOR_INT_MUL);
	ke_hash_add(sml, (fncp)&ke_vector_int_div, VECTOR_INT_DIV);
	ke_hash_add(sml, (fncp)&ke_vector_int_scale, VECTOR_INT_SCALE);
	ke_hash_add(sml, (fncp)&ke_vector_int_add_constant, VECTOR_INT_ADD_CONSTANT);
	ke_hash_add(sml, (fncp)&ke_vector_int_reverse, VECTOR_INT_REVERSE);
	ke_hash_add(sml, (fncp)&ke_vector_int_swap_elements, VECTOR_INT_SWAP_ELEMENTS);
	ke_hash_add(sml, (fncp)&ke_vector_int_memcpy, VECTOR_INT_MEMCPY);
	ke_hash_add(sml, (fncp)&ke_vector_int_swap, VECTOR_INT_SWAP);
	ke_hash_add(sml, (fncp)&ke_vector_int_min, VECTOR_INT_MIN);
	ke_hash_add(sml, (fncp)&ke_vector_int_max, VECTOR_INT_MAX);
	ke_hash_add(sml, (fncp)&ke_vector_int_isnull, VECTOR_INT_ISNULL);
	ke_hash_add(sml, (fncp)&ke_vector_int_ispos, VECTOR_INT_ISPOS);
	ke_hash_add(sml, (fncp)&ke_vector_int_isneg, VECTOR_INT_ISNEG);
	ke_hash_add(sml, (fncp)&ke_vector_int_isnonneg, VECTOR_INT_ISNONNEG);
	ke_hash_add(sml, (fncp)&ke_vector_int_equal, VECTOR_INT_EQUAL);
	ke_hash_add(sml, (fncp)&ke_vector_int_fscanf, VECTOR_INT_FSCANF);
	ke_hash_add(sml, (fncp)&ke_vector_int_fscanf, VECTOR_INT_SCAN);
	ke_hash_add(sml, (fncp)&ke_vector_int_fprintf, VECTOR_INT_FPRINTF);
	ke_hash_add(sml, (fncp)&ke_vector_int_fprintf, VECTOR_INT_PRINT);
	ke_hash_add(sml, (fncp)&ke_vector_int_fread, VECTOR_INT_FREAD);
	ke_hash_add(sml, (fncp)&ke_vector_int_fread, VECTOR_INT_READ);
	ke_hash_add(sml, (fncp)&ke_vector_int_fwrite, VECTOR_INT_FWRITE);
	ke_hash_add(sml, (fncp)&ke_vector_int_fwrite, VECTOR_INT_WRITE);

}

void ke_vector_print(sml_t* sml, token_t *k) {
    printf("Vector: %s\n", k->name);
    for(size_t i = 0; i < k->obj.vector->size; i++) {
        printf("%d : v:%g\n", (int)i, gsl_vector_get(k->obj.vector, i));
    }
}

void ke_vector_freemem(sml_t* sml,token_t *e) {
    if (e->obj.vector && e->vtype == KEV_VEC) {
        gsl_vector_free(e->obj.vector); ke_dec_memory(sml);
        e->obj.vector = NULL;
    }
}

void ke_vector_int_print(sml_t* sml,token_t *k) {
	printf("Vector: %s\n", k->name);
	for (size_t i = 0; i < k->obj.vector_int->size; i++) {
		printf("%d : v:%d\n", (int)i, gsl_vector_int_get(k->obj.vector_int, i));
	}
}

void ke_vector_int_freemem(sml_t* sml,token_t *e) {
	if (e->obj.vector_int && e->vtype == KEV_VEC_INT) {
		gsl_vector_int_free(e->obj.vector_int); ke_dec_memory(sml);
		e->obj.vector_int = NULL;
	}
}

max_stars_repo_name: vinej/sml

max_issues_repo_name: vinej/sml

max_forks_repo_name: vinej/sml


text:
// Copyright 2018, Kevin Lang, Oath Research

/*
gcc -Wall -pedantic -Iinclude -Llib -o multi_reading_linsolve multi_reading_linsolve.c -lgsl -lgslcblas
*/

#include <stdio.h>
#include <assert.h>
#include <gsl/gsl_linalg.h>

/*****************************************************************************/

#ifdef EXAMPLE_FOR_INTUITION

#define NUMEQNS 5
#define NUMVARS 4
double a_data[] =
    { 0.18, 0.60, 0.57, 0.96,
      0.41, 0.24, 0.99, 0.58,
      0.14, 0.30, 0.97, 0.66,
      0.51, 0.13, 0.19, 0.85,
      0.65, 0.43, 1.16, 1.51
    };

double b_data[] =
  { 1.0, 2.0, 3.0, 4.0,
    7.0
 };

#endif

/*****************************************************************************/

void read_a_data (int num_eqns, int num_vars, double * a_data) {
  int i, j, nxt;
  double tmp;
  nxt = 0;
  for (i = 0; i < num_eqns; i++) {
    for (j = 0; j < num_vars; j++) {
      int got = scanf ("%lf", &tmp);
      assert (got == 1);
      a_data[nxt++] = tmp;
    }
  }
  assert (nxt == num_eqns * num_vars);
}

/*****************************************************************************/

void read_b_data (int num_eqns, double * b_data) {
  int i, nxt;
  double tmp;
  nxt = 0;
  for (i = 0; i < num_eqns; i++) {
    int got = scanf ("%lf", &tmp);
    assert (got == 1);
    b_data[nxt++] = tmp;
  }
  assert (nxt == num_eqns);
}

/*****************************************************************************/

int main (int argc, char ** argv)
{
  if (argc != 3) {
    fprintf (stderr, "Usage: %s num_eqns num_vars < data_file\n", argv[0]);
    fprintf (stderr, "    The data_file contains a matrix followed by a target vector.\n");
    fprintf (stderr, "    The matrix should be in equation major order.\n");
    fprintf (stderr, "    This version solves multiple problems with matching matrix dimensions.\n");
    fprintf (stderr, "    Each problem should be preceding by a '1'.\n");
    fprintf (stderr, "    The last problem should be followed by a '0'.\n");
    return (-1);
  }

  int num_eqns  = atoi(argv[1]);
  int num_vars  = atoi(argv[2]);

  double * a_data = (double *) malloc (num_eqns * num_vars * sizeof(double));
  assert (a_data != NULL);

  double * b_data = (double *) malloc (num_eqns * sizeof(double));
  assert (b_data != NULL);

  while (1) {
    int keep_going = 0;
    int got = scanf ("%d", &keep_going);
    assert (got == 1);
    if (keep_going == 0) return 0;

    read_a_data (num_eqns, num_vars, a_data);
    read_b_data (num_eqns, b_data);

    gsl_matrix_view m = gsl_matrix_view_array (a_data, num_eqns, num_vars);
    gsl_vector *x     = gsl_vector_alloc (num_vars);
    gsl_vector *tau   = gsl_vector_alloc (num_vars);
    gsl_vector *residual = gsl_vector_alloc (num_eqns);
    gsl_vector_view b = gsl_vector_view_array (b_data, num_eqns);
    int s;

    s = gsl_linalg_QR_decomp (&m.matrix, tau);
    s = gsl_linalg_QR_lssolve (&m.matrix, tau, &b.vector, x, residual);
    gsl_vector_fprintf (stdout, x, "%.15g"); fflush (stdout);
    gsl_vector_free (x);
    gsl_vector_free (tau);
    gsl_vector_free (residual);
  }

  return 0;
}

    //  printf ("residual = \n"); gsl_vector_fprintf (stdout, residual, "%g");
    //      printf("solver got A: %d %d %.19g\n", i, j, tmp);
    //    printf("solver got B: %d %.19g\n", i, tmp);

ext: c

max_stars_repo_name: kevinjlang/cpc

max_issues_repo_name: kevinjlang/cpc

max_forks_repo_name: kevinjlang/cpc


text:
//STARTWHOLE
static char help[] = "Solve a 4x4 linear system using KSP.\n";

#include <petsc.h>

int main(int argc,char **args) {
    PetscErrorCode ierr;
    Vec        x, b;
    Mat        A;
    KSP        ksp;
    PetscInt   i, j[4] = {0, 1, 2, 3};                // j = column index
    PetscReal  ab[4] = {7.0, 1.0, 1.0, 3.0},          // vector entries
               aA[4][4] = {{ 1.0,  2.0,  3.0,  0.0},  // matrix entries
                           { 2.0,  1.0, -2.0, -3.0},
                           {-1.0,  1.0,  1.0,  0.0},
                           { 0.0,  1.0,  1.0, -1.0}};

    ierr = PetscInitialize(&argc,&args,NULL,help); if (ierr) return ierr;

    ierr = VecCreate(PETSC_COMM_WORLD,&b); CHKERRQ(ierr);
    ierr = VecSetSizes(b,PETSC_DECIDE,4); CHKERRQ(ierr);
    ierr = VecSetFromOptions(b); CHKERRQ(ierr);
    ierr = VecSetValues(b,4,j,ab,INSERT_VALUES); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(b); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(b); CHKERRQ(ierr);

    ierr = MatCreate(PETSC_COMM_WORLD,&A); CHKERRQ(ierr);
    ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,4,4); CHKERRQ(ierr);
    ierr = MatSetFromOptions(A); CHKERRQ(ierr);
    ierr = MatSetUp(A); CHKERRQ(ierr);
    for (i=0; i<4; i++) {   // set entries one row at a time
        ierr = MatSetValues(A,1,&i,4,j,aA[i],INSERT_VALUES); CHKERRQ(ierr);
    }
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);

    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,A,A); CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
    ierr = VecDuplicate(b,&x); CHKERRQ(ierr);
    ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);
    ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);

    KSPDestroy(&ksp);  MatDestroy(&A);
    VecDestroy(&x);  VecDestroy(&b);
    return PetscFinalize();
}
//ENDWHOLE


text:
#ifdef USE_MKL
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
#endif

// contains all calls to C-interface LAPACK functions (LAPACKE)
// https://www.netlib.org/lapack/lapacke.html

// returns the solution, x, to a real system of linear equations
//   A * x = b,
// solution is returned in b, i.e. b --> x, for return value 0
int specex_posv(int n, const double *A, const double *b){
  return LAPACKE_dposv(LAPACK_COL_MAJOR, 'L', n, 1, A, n, b, n);
}

// invert matrix A in place; A := inv(A)
// http://www.netlib.org/lapack/explore-html/d8/d63/dpotri_8f_source.html
int specex_potri(int n, const double *A){
  return LAPACKE_dpotri(LAPACK_COL_MAJOR, 'L', n, A, n);
}

max_stars_repo_name: tskisner/specex

max_issues_repo_name: tskisner/specex

max_forks_repo_name: tskisner/specex


text:
static char help[] =
"Time-dependent advection equation, in flux-conservative form, in 2D.\n"
"Option prefix -adv_.  Domain is (-1,1) x (-1,1).  Equation is\n"
"  u_t + div(a(x,y) u) = g(x,y,u).\n"
"Boundary conditions are periodic in x and y.  Cells are grid-point centered.\n"
"Allows flux-limited (non-oscillatory) method-of-lines discretization:\n"
"  none       O(h^1)  first-order upwinding (limiter = 0)\n"
"  centered   O(h^2)  linear centered\n"
"  vanleer    O(h^2)  van Leer (1974) limiter\n"
"  koren      O(h^3)  Koren (1993) limiter [default].\n"
"(There is separate control over the limiter in the residual and in the\n"
"Jacobian.  Only none and centered are implemented for the Jacobian.)\n"
"Solves either of two problems with initial conditions:\n"
"  straight   Figure 6.2, page 303, in Hundsdorfer & Verwer (2003) [default]\n"
"  rotation   Figure 20.5, page 461, in LeVeque (2002).\n"
"For straight, if final time is an integer and velocities are kept at default\n"
"values, then exact solution is known and L1,L2 errors are reported.\n\n";

#include <petsc.h>

//STARTCTX
typedef enum {STUMP, SMOOTH, CONE, BOX} InitialType;
static const char *InitialTypes[] = {"stump", "smooth", "cone", "box",
                                     "InitialType", "", NULL};

typedef enum {NONE, CENTERED, VANLEER, KOREN} LimiterType;
static const char *LimiterTypes[] = {"none","centered","vanleer","koren",
                                     "LimiterType", "", NULL};

typedef enum {STRAIGHT, ROTATION} ProblemType;
static const char *ProblemTypes[] = {"straight","rotation",
                                     "ProblemType", "", NULL};

typedef struct {
    ProblemType    problem;
    double         windx, windy,            // x,y velocity in STRAIGHT
                   (*initial_fcn)(double,double), // for STRAIGHT
                   (*limiter_fcn)(double),  // limiter used in RHS
                   (*jac_limiter_fcn)(double); // used in Jacobian
} AdvectCtx;
//ENDCTX

//STARTINITIAL
// equal to 1 in a disc of radius 0.2 around (-0.6,-0.6)
static double stump(double x, double y) {
    const double r = PetscSqrtReal((x+0.6)*(x+0.6) + (y+0.6)*(y+0.6));
    return (r < 0.2) ? 1.0 : 0.0;
}

// smooth (C^6) version of stump
static double smooth(double x, double y) {
    const double r = PetscSqrtReal((x+0.6)*(x+0.6) + (y+0.6)*(y+0.6));
    if (r < 0.2)
        return PetscPowReal(1.0 - PetscPowReal(r / 0.2,6.0),6.0);
    else
        return 0.0;
}

// cone of height 1 of base radius 0.35 centered at (-0.45,0.0)
static double cone(double x, double y) {
    const double r = PetscSqrtReal((x+0.45)*(x+0.45) + y*y);
    return (r < 0.35) ? 1.0 - r / 0.35 : 0.0;
}

// equal to 1 in square of side-length 0.5 (0.1,0.6) x (-0.25,0.25)
static double box(double x, double y) {
    if ((0.1 < x) && (x < 0.6) && (-0.25 < y) && (y < 0.25))
        return 1.0;
    else
        return 0.0;
}

static void* initialptr[] = {&stump, &smooth, &cone, &box};
//ENDINITIAL

//STARTLIMITERS
/* the centered-space method is linear */
static double centered(double theta) {
    return 0.5;
}

/* van Leer (1974) limiter is formula (1.11) in section III.1 of
Hundsdorfer & Verwer (2003) */
static double vanleer(double theta) {
    const double abstheta = PetscAbsReal(theta);
    return 0.5 * (theta + abstheta) / (1.0 + abstheta);
}

/* Koren (1993) limiter is formula (1.7) in section III.1 of
Hundsdorfer & Verwer (2003) */
static double koren(double theta) {
    const double z = (1.0/3.0) + (1.0/6.0) * theta;
    return PetscMax(0.0, PetscMin(1.0, PetscMin(z, theta)));
}

static void* limiterptr[] = {NULL, &centered, &vanleer, &koren};
//ENDLIMITERS

// velocity  a(x,y) = ( a^x(x,y), a^y(x,y) )
static double a_wind(double x, double y, int dir, AdvectCtx* user) {
    switch (user->problem) {
        case STRAIGHT:
            return (dir == 0) ? user->windx : user->windy;
        case ROTATION:
            return (dir == 0) ? y : - x;
        default:
            return 0.0;
    }
}

// source  g(x,y,u)
static double g_source(double x, double y, double u, AdvectCtx* user) {
    return 0.0;
}

//         d g(x,y,u) / d u
static double dg_source(double x, double y, double u, AdvectCtx* user) {
    return 0.0;
}

extern PetscErrorCode FormInitial(DMDALocalInfo*, Vec, AdvectCtx*);
extern PetscErrorCode DumpBinary(const char*, const char*, Vec);
extern PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo*, double,
        double**, double**, AdvectCtx*);
extern PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo*, double,
        double**, Mat, Mat, AdvectCtx*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    TS               ts;
    DM               da;
    Vec              u;
    DMDALocalInfo    info;
    double           hx, hy, t0, c, dt, tf;
    char             fileroot[PETSC_MAX_PATH_LEN] = "";
    int              steps;
    PetscBool        oneline = PETSC_FALSE, snesfdset, snesfdcolorset;
    InitialType      initial = STUMP;
    LimiterType      limiter = KOREN, jac_limiter = NONE;
    AdvectCtx        user;

    PetscInitialize(&argc,&argv,(char*)0,help);

    user.problem = STRAIGHT;
    user.windx = 2.0;
    user.windy = 2.0;
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
           "adv_", "options for advect.c", ""); CHKERRQ(ierr);
    ierr = PetscOptionsString("-dumpto","filename root for binary files with initial/final state",
           "advect.c",fileroot,fileroot,PETSC_MAX_PATH_LEN,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-initial",
           "shape of initial condition if problem==straight",
           "advect.c",InitialTypes,
           (PetscEnum)initial,(PetscEnum*)&initial,NULL); CHKERRQ(ierr);
    user.initial_fcn = initialptr[initial];
    ierr = PetscOptionsEnum("-limiter",
           "flux-limiter type used in RHS evaluation",
           "advect.c",LimiterTypes,
           (PetscEnum)limiter,(PetscEnum*)&limiter,NULL); CHKERRQ(ierr);
    user.limiter_fcn = limiterptr[limiter];
    ierr = PetscOptionsEnum("-jac_limiter",
           "flux-limiter type used in Jacobian (of RHS) evaluation",
           "advect.c",LimiterTypes,
           (PetscEnum)jac_limiter,(PetscEnum*)&jac_limiter,NULL); CHKERRQ(ierr);
    user.jac_limiter_fcn = limiterptr[jac_limiter];
    ierr = PetscOptionsEnum("-problem",
           "problem type",
           "advect.c",ProblemTypes,
           (PetscEnum)user.problem,(PetscEnum*)&user.problem,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-oneline",
           "in exact solution cases, show one-line output",
           "advect.c",oneline,&oneline,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-windx",
           "x component of wind for problem==straight",
           "advect.c",user.windx,&user.windx,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-windy",
           "y component of wind for problem==straight",
           "advect.c",user.windy,&user.windy,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    ierr = PetscOptionsHasName(NULL,NULL,"-snes_fd",&snesfdset); CHKERRQ(ierr);
    ierr = PetscOptionsHasName(NULL,NULL,"-snes_fd_color",&snesfdcolorset); CHKERRQ(ierr);
    if (snesfdset || snesfdcolorset) {
        user.jac_limiter_fcn = NULL;
        jac_limiter = 5;   // corresponds to empty string
    }

    ierr = DMDACreate2d(PETSC_COMM_WORLD,
               DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC,
               DMDA_STENCIL_STAR,              // no diagonal differencing
               5,5,PETSC_DECIDE,PETSC_DECIDE,  // default to hx=hx=0.2 grid
                                               //   (mx=my=5 allows -snes_fd_color)
               1, 2,                           // d.o.f & stencil width
               NULL,NULL,&da); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    hx = 2.0 / info.mx;  hy = 2.0 / info.my;
    ierr = DMDASetUniformCoordinates(da,    // grid is cell-centered
        -1.0+hx/2.0,1.0-hx/2.0,-1.0+hy/2.0,1.0-hy/2.0,0.0,1.0);CHKERRQ(ierr);

    ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
    ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
    ierr = TSSetDM(ts,da); CHKERRQ(ierr);
    ierr = DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
           (DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user); CHKERRQ(ierr);
    ierr = DMDATSSetRHSJacobianLocal(da,
           (DMDATSRHSJacobianLocal)FormRHSJacobianLocal,&user); CHKERRQ(ierr);
    ierr = TSSetType(ts,TSRK); CHKERRQ(ierr);  // defaults to -ts_rk_type 3bs

    // time axis: use CFL number of 0.5 to set initial time step, but note
    //            most methods adapt anyway
    if (user.problem == STRAIGHT)
        c = PetscMax(PetscAbsReal(user.windx)/hx, PetscAbsReal(user.windy)/hy);
    else
        c = PetscMax(1.0/hx, 1.0/hy);
    dt = 0.5 / c;
    ierr = TSSetTime(ts,0.0); CHKERRQ(ierr);
    ierr = TSSetMaxTime(ts,0.6); CHKERRQ(ierr);
    ierr = TSSetTimeStep(ts,dt); CHKERRQ(ierr);
    ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
    ierr = TSSetFromOptions(ts);CHKERRQ(ierr);

    ierr = DMCreateGlobalVector(da,&u); CHKERRQ(ierr);
    ierr = FormInitial(&info,u,&user); CHKERRQ(ierr);
    ierr = DumpBinary(fileroot,"_initial",u); CHKERRQ(ierr);
    ierr = TSGetTime(ts,&t0); CHKERRQ(ierr);
    ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);

    if (!oneline) {
        ierr = PetscPrintf(PETSC_COMM_WORLD,
               "solving problem %s with %s initial state on %d x %d grid,\n"
               "    cells dx=%g x dy=%g, limiter = %s, and jac_limiter = %s ...\n",
               ProblemTypes[user.problem],InitialTypes[initial],info.mx,info.my,
               hx,hy,LimiterTypes[limiter],LimiterTypes[jac_limiter]); CHKERRQ(ierr);
    }

    ierr = TSSolve(ts,u); CHKERRQ(ierr);

    ierr = TSGetStepNumber(ts,&steps); CHKERRQ(ierr);
    ierr = TSGetTime(ts,&tf); CHKERRQ(ierr);
    ierr = DumpBinary(fileroot,"_final",u); CHKERRQ(ierr);

    if (!oneline) {
        ierr = PetscPrintf(PETSC_COMM_WORLD,
                "completed %d steps to time %g\n",steps,tf); CHKERRQ(ierr);
    }

    if ( (user.problem == STRAIGHT) && (PetscAbs(fmod(tf+0.5e-8,1.0)) <= 1.0e-8)
         && (fmod(user.windx,2.0) == 0.0) && (fmod(user.windy,2.0) == 0.0) ) {
        // exact solution is same as initial condition
        Vec    uexact;
        double norms[2];
        ierr = VecDuplicate(u,&uexact); CHKERRQ(ierr);
        ierr = FormInitial(&info,uexact,&user); CHKERRQ(ierr);
        ierr = VecAXPY(u,-1.0,uexact); CHKERRQ(ierr); // u <- u + (-1.0) uexact
        ierr = VecNorm(u,NORM_1_AND_2,norms); CHKERRQ(ierr);
        norms[0] *= hx * hy;
        norms[1] *= PetscSqrtReal(hx * hy);
        VecDestroy(&uexact);
        if (oneline) {
            ierr = PetscPrintf(PETSC_COMM_WORLD,
                "%s,%s,%s,%d,%d,%g,%g,%d,%g,%.4e,%.4e\n",
                ProblemTypes[user.problem],InitialTypes[initial],
                LimiterTypes[limiter],info.mx,info.my,hx,hy,steps,tf,
                norms[0],norms[1]); CHKERRQ(ierr);
        } else {
            ierr = PetscPrintf(PETSC_COMM_WORLD,
                "errors |u-uexact|_{1,h} = %.4e, |u-uexact|_{2,h} = %.4e\n",
                norms[0],norms[1]); CHKERRQ(ierr);
        }
    }

    VecDestroy(&u);  TSDestroy(&ts);  DMDestroy(&da);
    return PetscFinalize();
}

PetscErrorCode FormInitial(DMDALocalInfo *info, Vec u, AdvectCtx* user) {
    PetscErrorCode ierr;
    int          i, j;
    double       hx, hy, x, y, **au;

    ierr = VecSet(u,0.0); CHKERRQ(ierr);  // clear it first
    ierr = DMDAVecGetArray(info->da, u, &au); CHKERRQ(ierr);
    hx = 2.0 / info->mx;  hy = 2.0 / info->my;
    for (j=info->ys; j<info->ys+info->ym; j++) {
        y = -1.0 + (j+0.5) * hy;
        for (i=info->xs; i<info->xs+info->xm; i++) {
            x = -1.0 + (i+0.5) * hx;
            switch (user->problem) {
                case STRAIGHT:
                    au[j][i] = (*user->initial_fcn)(x,y);
                    break;
                case ROTATION:
                    au[j][i] = cone(x,y) + box(x,y);
                    break;
                default:
                    SETERRQ(PETSC_COMM_WORLD,1,"invalid user->problem\n");
            }
        }
    }
    ierr = DMDAVecRestoreArray(info->da, u, &au); CHKERRQ(ierr);
    return 0;
}

// dumps to file; does nothing if string root is empty or NULL
PetscErrorCode DumpBinary(const char* root, const char* append, Vec u) {
    PetscErrorCode ierr;
    if ((root) && (strlen(root) > 0)) {
        PetscViewer  viewer;
        char         filename[PETSC_MAX_PATH_LEN] = "";
        sprintf(filename,"%s%s.dat",root,append);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
            "writing PETSC binary file %s ...\n",filename); CHKERRQ(ierr);
        ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,
            FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
        ierr = VecView(u,viewer); CHKERRQ(ierr);
        ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
    }
    return 0;
}

/* method-of-lines discretization gives ODE system  u' = G(t,u)
so our finite volume scheme computes
    G_ij = - (fluxE - fluxW)/hx - (fluxN - fluxS)/hy + g(x,y,U_ij)
but only east (E) and north (N) fluxes are computed
*/
//STARTFUNCTION
PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info, double t,
        double **au, double **aG, AdvectCtx *user) {
    int         i, j, q, dj, di;
    double      hx, hy, halfx, halfy, x, y, a,
                u_up, u_dn, u_far, theta, flux;

    // clear G first
    for (j = info->ys; j < info->ys + info->ym; j++)
        for (i = info->xs; i < info->xs + info->xm; i++)
            aG[j][i] = 0.0;
    // fluxes on cell boundaries are traversed in E,N order with indices
    // q=0 for E and q=1 for N; cell center has coordinates (x,y)
    hx = 2.0 / info->mx;  hy = 2.0 / info->my;
    halfx = hx / 2.0;     halfy = hy / 2.0;
    for (j = info->ys-1; j < info->ys + info->ym; j++) { // note -1 start
        y = -1.0 + (j+0.5) * hy;
        for (i = info->xs-1; i < info->xs + info->xm; i++) { // -1 start
            x = -1.0 + (i+0.5) * hx;
            if ((i >= info->xs) && (j >= info->ys)) {
                aG[j][i] += g_source(x,y,au[j][i],user);
            }
            for (q = 0; q < 2; q++) {   // E (q=0) and N (q=1) bdry fluxes
                if (q == 0 && j < info->ys)  continue;
                if (q == 1 && i < info->xs)  continue;
                di = 1 - q;
                dj = q;
                a = a_wind(x + halfx*di,y + halfy*dj,q,user);
                // first-order flux
                u_up = (a >= 0.0) ? au[j][i] : au[j+dj][i+di];
                flux = a * u_up;
                // use flux-limiter
                if (user->limiter_fcn != NULL) {
                    // formulas (1.2),(1.3),(1.6); H&V pp 216--217
                    u_dn = (a >= 0.0) ? au[j+dj][i+di] : au[j][i];
                    if (u_dn != u_up) {
                        u_far = (a >= 0.0) ? au[j-dj][i-di]
                                           : au[j+2*dj][i+2*di];
                        theta = (u_up - u_far) / (u_dn - u_up);
                        flux += a * (*user->limiter_fcn)(theta)
                                  * (u_dn-u_up);
                    }
                }
                // update owned G_ij on both sides of computed flux
                if (q == 0) {
                    if (i >= info->xs)
                        aG[j][i]   -= flux / hx;
                    if (i+1 < info->xs + info->xm)
                        aG[j][i+1] += flux / hx;
                } else {
                    if (j >= info->ys)
                        aG[j][i]   -= flux / hy;
                    if (j+1 < info->ys + info->ym)
                        aG[j+1][i] += flux / hy;
                }
            }
        }
    }
    return 0;
}
//ENDFUNCTION

PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo *info, double t,
        double **au, Mat J, Mat P, AdvectCtx *user) {
    PetscErrorCode ierr;
    const int   dir[4] = { 0, 1, 0, 1},  // use x (0) or y (1) component
                xsh[4] = { 1, 0,-1, 0},  ysh[4]   = { 0, 1, 0,-1};
    int         i, j, l, nc;
    double      hx, hy, halfx, halfy, x, y, a, v[9];
    MatStencil  col[9],row;

    ierr = MatZeroEntries(P); CHKERRQ(ierr);
    hx = 2.0 / info->mx;  hy = 2.0 / info->my;
    halfx = hx / 2.0;     halfy = hy / 2.0;
    for (j = info->ys; j < info->ys+info->ym; j++) {
        y = -1.0 + (j+0.5) * hy;
        row.j = j;
        for (i = info->xs; i < info->xs+info->xm; i++) {
            x = -1.0 + (i+0.5) * hx;
            row.i = i;
            col[0].j = j;  col[0].i = i;
            v[0] = dg_source(x,y,au[j][i],user);
            nc = 1;
            for (l = 0; l < 4; l++) {   // loop over cell boundaries: E, N, W, S
                a = a_wind(x + halfx*xsh[l],y + halfy*ysh[l],dir[l],user);
                if (user->jac_limiter_fcn == NULL) {
                    // Jacobian is from upwind fluxes
                    switch (l) {
                        case 0:
                            col[nc].j = j;
                            col[nc].i = (a >= 0.0) ? i : i+1;
                            v[nc++] = - a / hx;
                            break;
                        case 1:
                            col[nc].j = (a >= 0.0) ? j : j+1;
                            col[nc].i = i;
                            v[nc++] = - a / hy;
                            break;
                        case 2:
                            col[nc].j = j;
                            col[nc].i = (a >= 0.0) ? i-1 : i;
                            v[nc++] = a / hx;
                            break;
                        case 3:
                            col[nc].j = (a >= 0.0) ? j-1 : j;
                            col[nc].i = i;
                            v[nc++] = a / hy;
                            break;
                    }
                } else if (user->jac_limiter_fcn == &centered) {
                    // Jacobian is from centered fluxes
                    switch (l) {
                        case 0:
                            col[nc].j = j;  col[nc].i = i;    v[nc++] = - a / (2.0*hx);
                            col[nc].j = j;  col[nc].i = i+1;  v[nc++] = - a / (2.0*hx);
                            break;
                        case 1:
                            col[nc].j = j;    col[nc].i = i;  v[nc++] = - a / (2.0*hy);
                            col[nc].j = j+1;  col[nc].i = i;  v[nc++] = - a / (2.0*hy);
                            break;
                        case 2:
                            col[nc].j = j;  col[nc].i = i-1;  v[nc++] = a / (2.0*hx);
                            col[nc].j = j;  col[nc].i = i;    v[nc++] = a / (2.0*hx);
                            break;
                        case 3:
                            col[nc].j = j-1;  col[nc].i = i;  v[nc++] = a / (2.0*hy);
                            col[nc].j = j;    col[nc].i = i;  v[nc++] = a / (2.0*hy);
                            break;
                    }
                } else {
                    SETERRQ(PETSC_COMM_WORLD,1,"only Jacobian cases none|centered are implemented\n");
                }
            }
            ierr = MatSetValuesStencil(P,1,&row,nc,col,v,ADD_VALUES); CHKERRQ(ierr);
        }
    }
    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}


text:
/* matrix/gsl_matrix_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_MATRIX_DOUBLE_H__
#define __GSL_MATRIX_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_double.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  double * data;
  gsl_block * block;
  int owner;
} gsl_matrix;

typedef struct
{
  gsl_matrix matrix;
} _gsl_matrix_view;

typedef _gsl_matrix_view gsl_matrix_view;

typedef struct
{
  gsl_matrix matrix;
} _gsl_matrix_const_view;

typedef const _gsl_matrix_const_view gsl_matrix_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix *
gsl_matrix_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix *
gsl_matrix_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix *
gsl_matrix_alloc_from_block (gsl_block * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_EXPORT
gsl_matrix *
gsl_matrix_alloc_from_matrix (gsl_matrix * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
gsl_vector *
gsl_vector_alloc_row_from_matrix (gsl_matrix * m,
                                        const size_t i);

GSL_EXPORT
gsl_vector *
gsl_vector_alloc_col_from_matrix (gsl_matrix * m,
                                        const size_t j);

GSL_EXPORT void gsl_matrix_free (gsl_matrix * m);

/* Views */

GSL_EXPORT
_gsl_matrix_view
gsl_matrix_submatrix (gsl_matrix * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_view
gsl_matrix_row (gsl_matrix * m, const size_t i);

GSL_EXPORT
_gsl_vector_view
gsl_matrix_column (gsl_matrix * m, const size_t j);

GSL_EXPORT
_gsl_vector_view
gsl_matrix_diagonal (gsl_matrix * m);

GSL_EXPORT
_gsl_vector_view
gsl_matrix_subdiagonal (gsl_matrix * m, const size_t k);

GSL_EXPORT
_gsl_vector_view
gsl_matrix_superdiagonal (gsl_matrix * m, const size_t k);

GSL_EXPORT
_gsl_matrix_view
gsl_matrix_view_array (double * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_view
gsl_matrix_view_array_with_tda (double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_view
gsl_matrix_view_vector (gsl_vector * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_view
gsl_matrix_view_vector_with_tda (gsl_vector * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_const_view
gsl_matrix_const_submatrix (const gsl_matrix * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_const_view
gsl_matrix_const_row (const gsl_matrix * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_const_view
gsl_matrix_const_column (const gsl_matrix * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_const_view
gsl_matrix_const_diagonal (const gsl_matrix * m);

GSL_EXPORT
_gsl_vector_const_view
gsl_matrix_const_subdiagonal (const gsl_matrix * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_const_view
gsl_matrix_const_superdiagonal (const gsl_matrix * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_const_view
gsl_matrix_const_view_array (const double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_const_view
gsl_matrix_const_view_array_with_tda (const double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_const_view
gsl_matrix_const_view_vector (const gsl_vector * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_const_view
gsl_matrix_const_view_vector_with_tda (const gsl_vector * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT double   gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j);
GSL_EXPORT void    gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x);

GSL_EXPORT double * gsl_matrix_ptr(gsl_matrix * m, const size_t i, const size_t j);
GSL_EXPORT const double * gsl_matrix_const_ptr(const gsl_matrix * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_set_zero (gsl_matrix * m);
GSL_EXPORT void gsl_matrix_set_identity (gsl_matrix * m);
GSL_EXPORT void gsl_matrix_set_all (gsl_matrix * m, double x);

GSL_EXPORT int gsl_matrix_fread (FILE * stream, gsl_matrix * m) ;
GSL_EXPORT int gsl_matrix_fwrite (FILE * stream, const gsl_matrix * m) ;
GSL_EXPORT int gsl_matrix_fscanf (FILE * stream, gsl_matrix * m);
GSL_EXPORT int gsl_matrix_fprintf (FILE * stream, const gsl_matrix * m, const char * format);

GSL_EXPORT int gsl_matrix_memcpy(gsl_matrix * dest, const gsl_matrix * src);
GSL_EXPORT int gsl_matrix_swap(gsl_matrix * m1, gsl_matrix * m2);

GSL_EXPORT int gsl_matrix_swap_rows(gsl_matrix * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_swap_columns(gsl_matrix * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_swap_rowcol(gsl_matrix * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_transpose (gsl_matrix * m);
GSL_EXPORT int gsl_matrix_transpose_memcpy (gsl_matrix * dest, const gsl_matrix * src);

GSL_EXPORT double gsl_matrix_max (const gsl_matrix * m);
GSL_EXPORT double gsl_matrix_min (const gsl_matrix * m);
GSL_EXPORT void gsl_matrix_minmax (const gsl_matrix * m, double * min_out, double * max_out);

GSL_EXPORT void gsl_matrix_max_index (const gsl_matrix * m, size_t * imax, size_t *jmax);
GSL_EXPORT void gsl_matrix_min_index (const gsl_matrix * m, size_t * imin, size_t *jmin);
GSL_EXPORT void gsl_matrix_minmax_index (const gsl_matrix * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_EXPORT int gsl_matrix_isnull (const gsl_matrix * m);

GSL_EXPORT int gsl_matrix_add (gsl_matrix * a, const gsl_matrix * b);
GSL_EXPORT int gsl_matrix_sub (gsl_matrix * a, const gsl_matrix * b);
GSL_EXPORT int gsl_matrix_mul_elements (gsl_matrix * a, const gsl_matrix * b);
GSL_EXPORT int gsl_matrix_div_elements (gsl_matrix * a, const gsl_matrix * b);
GSL_EXPORT int gsl_matrix_scale (gsl_matrix * a, const double x);
GSL_EXPORT int gsl_matrix_add_constant (gsl_matrix * a, const double x);
GSL_EXPORT int gsl_matrix_add_diagonal (gsl_matrix * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_get_row(gsl_vector * v, const gsl_matrix * m, const size_t i);
GSL_EXPORT int gsl_matrix_get_col(gsl_vector * v, const gsl_matrix * m, const size_t j);
GSL_EXPORT int gsl_matrix_set_row(gsl_matrix * m, const size_t i, const gsl_vector * v);
GSL_EXPORT int gsl_matrix_set_col(gsl_matrix * m, const size_t j, const gsl_vector * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
double
gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
double *
gsl_matrix_ptr(gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (double *) (m->data + (i * m->tda + j)) ;
}

extern inline
const double *
gsl_matrix_const_ptr(const gsl_matrix * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const double *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_DOUBLE_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
// cosmo.h
// written by Suman Bhattacharya (https://github.com/bsuman79), modified by Samuel Flender
// defines the cosmology and allows calculations of various cosmological params

#ifndef cosmo_Header_included
#define cosmo_Header_included

#include <stdio.h>
#include <math.h>
#include <iostream>
#include <stdlib.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>

//double angdiam_func(double x, void * p);

using namespace std;

class cosmo {

protected:
	float H0, Omega_M, rho_crit, a, Omega_k, Omega_L, wt;
	float PI, m_sun, G, mpc, gnewt;

public:

	cosmo(float inp1, float inp2, float inp4, float inp5)
	{
		// cosmological params at z =0
		H0 = inp1;
		Omega_M = inp2;
		Omega_k = inp4;
		Omega_L = 1.0- Omega_M- Omega_k;
		wt= inp5;
		PI = 4.*atan(1.);
		m_sun = 1.98892e30;
		G = 6.67e-11/1.0e9; // in km^3 kg^-1 s^-2
		mpc = 3.0857e19; // in km
		gnewt = G*m_sun/mpc; //for cosm params (in km^2 Mpc msun^-1 s^-2)
	}

	float get_H0() {
		return H0;
	}

	float get_hubble70() {
		return H0/70.0;
	}

	float scale_fact(float redshift) {
		a = 1.0/(1.0+redshift);
		return a;
	}

	float hubblez(float redshift) {
		scale_fact(redshift);
		return H0 * Efact(redshift);
	}

	float calc_rho_crit(float redshift) {
		float Hz = hubblez(redshift);
		rho_crit = pow(Hz,2)*3./(8.*PI*gnewt);
		return rho_crit;
	}

	float Omega_Mz(float redshift) {
		scale_fact(redshift);
		float Hz = hubblez(redshift);
		return (Omega_M/pow(a,3))*pow(H0/Hz,2);
	}

	float Delta_vir(float redshift) {
		float x;
		x  = Omega_Mz(redshift) - 1.0;
		return (18*pow(PI,2)) + (82*x) - (39*pow(x,2));
	}

	float Efact(float redshift) {
		scale_fact(redshift);
		return sqrt(Omega_M/pow(a,3) + Omega_L/pow(a, 3*wt+3) + (1.0-Omega_M-Omega_L)/pow(a,2));
	}

	float ang_diam(float redshift) {
		// note rcutoff is in units of Rvir
		int dummy;
		float result;
		double x[1000], y[1000];
		float speed_of_light = 2.9979e5; //km/s
		double units = scale_fact(redshift)*speed_of_light/H0;

		for(dummy=0;dummy<1000;dummy++){
			x[dummy]= 0.0+ redshift/999*dummy;
			y[dummy] = 1.0/(sqrt( Omega_M*pow(1+x[dummy],3) + Omega_L*pow(1+x[dummy], 3*wt+3) + (1.0-Omega_M-Omega_L)*pow(1+x[dummy],2)));
		}
		if (redshift==0.0) return 0.0;
		else {
			gsl_interp_accel *acc
			= gsl_interp_accel_alloc ();

			gsl_spline *spline
			= gsl_spline_alloc (gsl_interp_cspline, 1000);

			gsl_spline_init (spline, x, y, 1000);
			result=units*gsl_spline_eval_integ (spline, x[0], x[dummy-1], acc);
			gsl_spline_free (spline);
			gsl_interp_accel_free (acc);
		}
		return result;
	}

	float comov_dist(float redshift) { //comoving distance to redshift z in Mpc
		// note rcutoff is in units of Rvir
		int dummy;
		float result;
		double x[1000], y[1000];
		float speed_of_light = 2.9979e5; //km/s
		double units = speed_of_light/H0;

		for(dummy=0;dummy<1000;dummy++){
			x[dummy]= 0.0+ redshift/999*dummy;
			y[dummy] = 1.0/(sqrt( Omega_M*pow(1+x[dummy],3) + Omega_L*pow(1+x[dummy], 3*wt+3) + (1.0-Omega_M-Omega_L)*pow(1+x[dummy],2)));
		}
		if (redshift==0.0) return 0.0;
		else {
			gsl_interp_accel *acc
			= gsl_interp_accel_alloc ();

			gsl_spline *spline
			= gsl_spline_alloc (gsl_interp_cspline, 1000);

			gsl_spline_init (spline, x, y, 1000);
			result=units*gsl_spline_eval_integ (spline, x[0], x[dummy-1], acc);
			gsl_spline_free (spline);
			gsl_interp_accel_free (acc);
		}
		return result;
	}

	float cosmic_time(float redshift) {
		// function to calculate cosmic time to redshift given cosmological model
		//good to ~1% for z<10

		// these are weight functions for solving the differential equation
		float x[16] = {0.0950,0.2816,0.4580,0.6179,0.7554,0.8656,0.9446,0.9894,-0.0950,-0.2816,-0.4580,-0.6179,-0.7554,-0.8656,-0.9446,-0.9894};
		float w[16] = {0.1895,0.1826,0.1692,0.1496,0.1246,0.0952,0.0623,0.0272,0.1895,0.1826,0.1692,0.1496,0.1246,0.0952,0.0623,0.0272};

		float a0 = 0.0, a1, Delta_a, h2, sum, cosm_h0, den0 = Omega_M, den1 = Omega_L, den2;
		float ss, ak, adot, cosm_t;
		float time_units, secperyr = 3.15569260e7, age_of_univ  = 8.93120195E+17/2.0; // in seconds
		int Nsteps, i, k;
		Nsteps = 10000;
		a1 = 1.0/(1.0+redshift);
		Delta_a = (a1-a0)/Nsteps;
		h2 = Delta_a/2.0;
		sum = 0.0;
		cosm_h0 = sqrt(8.0*PI/3.0);

		time_units = (age_of_univ*2/secperyr/1.0E9)*(1.0E2/H0); // converts to Gyr
		// now solve friedman eqn.
		for (i=0;i<Nsteps;i++) {
			ss = 0.0;
			for (k=0;k<16;k++) {
				ak = a0+h2*(1.0+x[k]);
				den2 = den1*dynrho(ak,wt);
				adot = sqrt( den0/ak + den2*ak*ak + 1.0-den0-den1 );
				ss = ss + w[k]/adot;
			}
			ss *= h2;
			a0 += Delta_a;
			sum += ss;
		}
		cosm_t = sum/cosm_h0;
		cosm_t = cosm_t*time_units;
		return cosm_t;
	}

	float dynrho(float a, float wdyn) {
		//       CALDWELL: This function RETURNs the energy density
		//       in the dynamical phi field, given the scale factor
		//       such that at a=1, dynrho = 1
		int n_dyn_q1,n_dyn_q2,n_dyn_q3;
		float drho;
		if (wdyn==-1.0) {
			n_dyn_q1 = 0;
			n_dyn_q2 = 0;
			n_dyn_q3 = 0;
		}
		else {
			n_dyn_q1 = 1;
			n_dyn_q2 = 1;
			n_dyn_q3 = 0;
		}
		//  did I ask for dynamical field? Leave this as an option for later
		if (n_dyn_q1==1) {
			if(n_dyn_q2==1) drho = pow(a,(-3*(1.0 + wdyn)));
		}
		//  otherwise, set the energy density to 1
		else  drho = 1.0;
		return drho;
	}

	friend class cluster;
	friend class growth;
	//friend double angdiam_func(double x, void * p);
};

/*double angdiam_func(double x, void * p) {
 double *params = (double *) p;
 double OmegaM = params[0];
 double OmegaL = params[1];
 double wt= params[2];
 double y, a = 1.0/(1.0+x);
 y = 1.0/(sqrt( OmegaM/pow(a,3) + OmegaL/pow(a, 3*wt+3) + (1.0-OmegaM-OmegaL)/pow(a,2)));
 return y;
 }*/
#endif

max_stars_repo_name: sflender/pairwise

max_issues_repo_name: sflender/pairwise

max_forks_repo_name: sflender/pairwise


text:
/*
*
*  Copyright 2011, 2012 by the CALATK development team
*
*   Licensed under the Apache License, Version 2.0 (the "License");
*   you may not use this file except in compliance with the License.
*   You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
*   Unless required by applicable law or agreed to in writing, software
*   distributed under the License is distributed on an "AS IS" BASIS,
*   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
*   See the License for the specific language governing permissions and
*   limitations under the License.
*
*
*/

#ifndef C_SOLVER_NLOPT_H
#define C_SOLVER_NLOPT_H

#include <string>
#include <math.h>
#include <nlopt.h>

namespace CALATK
{

/**
  * Implements an LBFGS solver
  *
  */
template < class TState >
class CSolverNLOpt
    : public CSolver< TState >
{
public:
  typedef typename TState::TFloat T;
  typedef CSolver< TState > Superclass;
  typedef typename Superclass::ObjectiveFunctionType ObjectiveFunctionType;
  typedef typename Superclass::CEnergyValues CEnergyValues;

  CSolverNLOpt();
  ~CSolverNLOpt();

  bool SolvePreInitialized();

  // specializations
  bool SolvePreInitialized( double* ptr, long int liNumberOfStateVectorElements );
  bool SolvePreInitialized( float* ptr, long int liNumberOfStateVectorElements );

  virtual void SetAutoConfiguration( CJSONConfiguration * combined, CJSONConfiguration * cleaned );

  SetMacro( initial_step1, T );
  GetMacro( initial_step1, T );

  SetMacro( maxeval, unsigned int );
  GetMacro( maxeval, unsigned int );

  SetMacro( xtol_rel, T );
  GetMacro( xtol_rel, T );

  SetMacro( ftol_rel, T );
  GetMacro( ftol_rel, T );

  SetMacro( vector_storage, unsigned int );
  GetMacro( vector_storage, unsigned int );

protected:

  struct SUserData
  {
    CSolverNLOpt< TState > *pSolver;
  };

  std::string GetStringForStatusCode( int ret );

  T* m_x;

    // specialization for double
   static double _evaluate(
       unsigned n,
       const double *x,
       double *g,
       void *my_func_data
       );

    // specialization for float
    static float _evaluate(
        unsigned n,
        const float *x,
        float *g,
        void *my_func_data
        );

     // generic evaluate and progress methods

    T evaluate(
        unsigned n,
        const T *x,
        T *g,
        void *my_func_data
        );

    T m_initial_step1;
    const T Defaultinitial_step1;
    bool m_ExternallySetinitial_step1;

    unsigned int m_maxeval;
    const unsigned int Defaultmaxeval;
    bool m_ExternallySetmaxeval;

    T m_xtol_rel;
    const T Defaultxtol_rel;
    bool m_ExternallySetxtol_rel;

    T m_ftol_rel;
    const T Defaultftol_rel;
    bool m_ExternallySetftol_rel;

    unsigned int m_vector_storage;
    const unsigned int Defaultvector_storage;
    bool m_ExternallySetvector_storage;

};

#include "CSolverNLOpt.txx"

} // end namespace

#endif // C_SOLVER_NLOPT_H

max_issues_repo_name: cpatrick/calatk

max_forks_repo_name: cpatrick/calatk


text:
#ifndef STK_Matrix_h
#define STK_Matrix_h

#include "common.h"
#include "Error.h"
#include "BasicVector.h"

#include <stddef.h>
#include <stdlib.h>
#include <stdexcept>
#include <iostream>

#ifdef HAVE_ATLAS
extern "C"{
  #ifdef USE_MKL
    #include "mkl.h"
  #else
    #include <cblas.h>
    #include <atlas/clapack.h>
  #endif
}
#endif

//#define TRACE_MATRIX_OPERATIONS
#define CHECKSIZE

namespace STK
{

  //  class matrix_error : public std::logic_error {};/
  //  class matrix_sizes_error : public matrix_error {};

  /// defines a storage type
  typedef enum
  {
    STORAGE_UNDEFINED = 0,
    STORAGE_REGULAR,
    STORAGE_TRANSPOSED
  } StorageType;

  // declare the class so the header knows about it
  template<typename _ElemT> class Matrix;
  template<typename _ElemT> class MatrixRange;

  // we need to declare the friend << operator here
  template<typename _ElemT>
    std::ostream & operator << (std::ostream & out, Matrix<_ElemT> & m);

  /** **************************************************************************
   ** **************************************************************************
   *  @brief Provides a matrix class
   *
   *  This class provides a way to work with matrices in STK.
   *  It encapsulates basic operations and memory optimizations.
   *
   */
  template<typename _ElemT>
    class Matrix
    {
    public:
      // HTK parameter file header (see HTK manual)
      typedef struct
      {
        int   nSamples;
        int   sampPeriod;
        short sampSize;
        short paramKind;
      } HtkHeader;

      /// defines a type of this
      typedef Matrix<_ElemT>    ThisType;

      // Constructors

      /// Empty constructor
      Matrix<_ElemT> ():
        mMRows(0), mMCols(0), mStride(0), mpData(NULL)
#ifdef STK_MEMALIGN_MANUAL
        , mpFreeData(NULL)
#endif
      {}

      /// Copy constructor
      Matrix<_ElemT> (const ThisType & t);

      /// Basic constructor
      Matrix<_ElemT> (const size_t r,
                      const size_t c);

      /// Destructor
      virtual
      ~Matrix<_ElemT> ();

      /// Initializes matrix (if not done by constructor)
      void
      Init(const size_t r,
           const size_t c);

      /**
       * @brief Dealocates the matrix from memory and resets the dimensions to (0, 0)
       */
      void
      Destroy();

      /**
       * @brief Returns @c true if matrix is initialized
       */
      const bool
      IsInitialized() const
      { return mpData != NULL; }

      /// Returns number of rows in the matrix
      const size_t
      Rows() const
      {
        return mMRows;
      }

      /// Returns number of columns in the matrix
      const size_t
      Cols() const
      {
        return mMCols;
      }

      /// Returns number of columns in the matrix memory
      const size_t
      Stride() const
      {
        return mStride;
      }

      /**
       *  @brief Gives access to a specified matrix row without range check
       *  @return Pointer to the const array
       */
      const _ElemT*
      cpData () const
      {
        return mpData;
      }

      /**
       *  @brief Gives access to a specified matrix row without range check
       *  @return Pointer to the non-const data array
       */
      _ElemT* __attribute__((aligned(16)))
      pData () const
      {
        return mpData;
      }

      /// Returns size of matrix in memory
      const size_t
      MSize() const
      {
        return mMRows * mStride * sizeof(_ElemT);
      }

      bool
      LoadHTK(const char* pFileName);

      //########################################################################
      //########################################################################
      // Math stuff
      /**
       * @brief Performs diagonal scaling
       * @param pDiagVector Array representing matrix diagonal
       * @return Refference to this
       */
      ThisType&
      DiagScale(const _ElemT* pDiagVector);

      ThisType&
      DiagScale(const BasicVector<_ElemT>& rDiagVector);

      /**
       *  @brief Performs vector multiplication on a and b and and adds the
       *         result to this (elem by elem)
       */
      ThisType &
      AddCMtMMul(const _ElemT& c, const ThisType& a, const ThisType& b);

      ThisType &
      AddMMMul(const ThisType& a, const ThisType& b);

      ThisType &
      AddMMtMul(const ThisType& a, const ThisType& b);

      ThisType &
      AddCMMtMul(const _ElemT& c, const ThisType& a, const ThisType& b);

      ThisType &
      AddCMMul(const _ElemT& c, const ThisType& a);

      ThisType &
      RepMMSub(const ThisType& a, const ThisType& b);

      ThisType &
      RepMMMul(const ThisType& a, const ThisType& b);

      ThisType &
      RepMtMMul(const ThisType& a, const ThisType& b);

      ThisType &
      AddCVVtMul(const _ElemT& c, const BasicVector<_ElemT>& rA,
                 const BasicVector<_ElemT>& rB);

      ThisType &
      AddCVVtMul(const _ElemT& c, const _ElemT* pA, const _ElemT* pB);

      ThisType &
      AddCVVt(const _ElemT& c, const BasicVector<_ElemT>& rA, const _ElemT* pB);

      ThisType &
      DivC(const _ElemT& c);

      //########################################################################
      //########################################################################

      ThisType &
      DotMul(const ThisType& a);

      ThisType&
      NormalizeRows();

      /**
       *  @brief Returns sum of all elements
       */
      const _ElemT&
      Sum() const;

      /**
       *  @brief Performs log on all elements
       *         result to this (elem by elem)
       */
      ThisType&
      Log();

      /**
       *  @brief Performs fast sigmoid on row vectors
       *         result to this (elem by elem)
       */
      ThisType &
      FastRowSigmoid();

      /**
       *  @brief Performs fast softmax on row vectors
       *         result to this (elem by elem)
       */
      ThisType &
      FastRowSoftmax();

      /**
       *  @brief Performs sigmoid on row vectors
       *         result to this (elem by elem)
       */
      ThisType &
      RowSigmoid();

      /**
       *  @brief Performs softmax on row vectors
       *         result to this (elem by elem)
       */
      ThisType &
      RowSoftmax();

      /**
       *  @brief Performs matrix inversion
       */
      ThisType &
      Invert();

      ThisType &
      Clear();

      /**
       *  @brief Gives access to a specified matrix row without range check
       *  @return pointer to the first field of the row
       */
      inline _ElemT*
      operator []  (size_t i) const
      {
        return mpData + (i * mStride);
      }

      /**
       *  @brief Gives access to a specified matrix row with range check
       *  @return pointer to the first field of the row
       */
      _ElemT*
      Row (const size_t r)
      {
        if (0 <= r && r < mMRows)
        {
          return this->operator[] (r);
        }
        else
        {
          throw std::out_of_range("Matrix row out of range");
        }
      }

      /**
       *  @brief Gives access to matrix elements (row, col)
       *  @return pointer to the desired field
       */
      _ElemT&
      operator () (const size_t r, const size_t c)
      { return *(mpData + r * mStride + c); }

      friend std::ostream &
      operator << <> (std::ostream & out, ThisType & m);

      void PrintOut(char *file);
      void ReadIn(char *file);

      /**
       * @brief Returns a matrix sub-range
       * @param ro Row offset
       * @param r  Rows in range
       * @param co Column offset
       * @param c  Coluns in range
       * See @c MatrixRange class for details
       */
      MatrixRange<_ElemT>
      Range(const size_t    ro, const size_t    r,
            const size_t    co, const size_t    c)
      { return MatrixRange<_ElemT>(*this, ro, r, co, c); }

    protected:
      inline void swap4b(void *a);
      inline void swap2b(void *a);

    protected:

      /// these atributes store the real matrix size as it is stored in memory
      /// including memalignment
      size_t    mMRows;       ///< Number of rows
      size_t    mMCols;       ///< Number of columns
      size_t    mStride;      ///< true number of columns for the internal matrix.
                              ///< This number may differ from M_cols as memory
                              ///< alignment might be used
      /// data memory area
      _ElemT*   mpData;

#ifdef STK_MEMALIGN_MANUAL
      /// data to be freed (in case of manual memalignment use, see common.h)
      _ElemT*   mpFreeData;
#endif
    }; // class Matrix

  /** **************************************************************************
   ** **************************************************************************
   *  @brief Provides a window matrix abstraction class
   *
   *  This class provides a way to work with matrix cutouts in STK.
   *  It encapsulates basic operations and memory optimizations.
   *
   */
  template<typename _ElemT>
    class WindowMatrix : public Matrix<_ElemT>
    {
    protected:
      /// points to the original begining of the data array
      /// The data atribute points now to the begining of the window
      _ElemT * mpOrigData;

      //@{
      /// these atributes store the real matrix size as it is stored in memory
      /// including memalignment
      size_t  mOrigMRows;       ///< Number of rows
      size_t  mOrigMCols;       ///< Number of columns
      size_t  mOrigMRealCols;   ///< true number of columns for the internal matrix.
                                ///< This number may differ from M_cols as memory
                                ///< alignment might be used
      size_t  mTRowOff;         ///< First row of the window
      size_t  mTColOff;         ///< First column of the window
      //@}

    public:
      /// defines a type of this
      typedef WindowMatrix<_ElemT>    ThisType;

      /// Empty constructor
      WindowMatrix() : Matrix<_ElemT>() {};

      /// Copy constructor
      WindowMatrix(const ThisType & rT);

      /// Basic constructor
      WindowMatrix(const size_t r, const size_t c):
        Matrix<_ElemT>(r, c), // create the base class
        mOrigMRows    (Matrix<_ElemT>::mMRows),
        mOrigMCols    (Matrix<_ElemT>::mMCols),
        mOrigMRealCols(Matrix<_ElemT>::mStride),
        mTRowOff(0), mTColOff(0)          // set the offset
      {
        mpOrigData = Matrix<_ElemT>::mpData;
      }

      /// The destructor
      ~WindowMatrix<_ElemT>()
      {
        Matrix<_ElemT>::mpData = this->mpOrigData;
      }

      /// sets the size of the window (whole rows)
      void
      SetSize(const size_t ro,
              const size_t r);

      /// sets the size of the window
      void
      SetSize(const size_t ro,
              const size_t r,
              const size_t co,
              const size_t c);

      /**
       *  @brief Resets the window to the default (full) size
       */
      void
      Reset ()
      {
        Matrix<_ElemT>::mpData    =  mpOrigData;
        Matrix<_ElemT>::mMRows    =  mOrigMRows;
        Matrix<_ElemT>::mMCols    =  mOrigMCols;
        Matrix<_ElemT>::mStride   =  mOrigMRealCols;

        mTRowOff = 0;
        mTColOff = 0;
      }
    };

  /** **************************************************************************
   ** **************************************************************************
   *  @brief Sub-matrix representation
   *
   *  This class provides a way to work with matrix cutouts in STK.
   *
   *
   */
  template<typename _ElemT>
    class MatrixRange : public Matrix<_ElemT>
    {
    public:
      /// Constructor
      MatrixRange(const Matrix<_ElemT>& rT,
                  const size_t    ro,
                  const size_t    r,
                  const size_t    co,
                  const size_t    c);

      /// The destructor
      virtual
      ~MatrixRange<_ElemT>()
      {
#ifndef STK_MEMALIGN_MANUAL
        Matrix<_ElemT>::mpData = NULL;
#else
        Matrix<_ElemT>::mpFreeData = NULL;
#endif
      }
    };
} // namespace STK

//*****************************************************************************
//*****************************************************************************
// we need to include the implementation
#include "Matrix.tcc"
//*****************************************************************************
//*****************************************************************************

/******************************************************************************
 ******************************************************************************
 * The following section contains specialized template definitions
 * whose implementation is in Matrix.cc
 */

namespace STK
{
  template<>
    Matrix<float> &
    Matrix<float>::
    FastRowSigmoid();

  template<>
    Matrix<double> &
    Matrix<double>::
    FastRowSigmoid();

  template<>
    Matrix<float> &
    Matrix<float>::
    RowSigmoid();

  template<>
    Matrix<double> &
    Matrix<double>::
    RowSigmoid();

  template<>
    Matrix<float> &
    Matrix<float>::
    DiagScale(const BasicVector<float>& rDiagVector);

  template<>
    Matrix<float> &
    Matrix<float>::
    DiagScale(const float* pDiagVector);

  template<>
    Matrix<float> &
    Matrix<float>::
    FastRowSoftmax();

  template<>
    Matrix<double> &
    Matrix<double>::
    FastRowSoftmax();

  template<>
    Matrix<float> &
    Matrix<float>::
    RowSoftmax();

  template<>
    Matrix<double> &
    Matrix<double>::
    RowSoftmax();

  template<>
    Matrix<float> &
    Matrix<float>::
    AddCMtMMul(const float& c, const Matrix<float>& a,
                               const Matrix<float>& b);

  template<>
    Matrix<double> &
    Matrix<double>::
    AddCMtMMul(const double& c, const Matrix<double>& a,
                                const Matrix<double>& b);

  template<>
    Matrix<float> &
    Matrix<float>::
    AddMMMul(const Matrix<float>& a, const Matrix<float>& b);

  template<>
    Matrix<double> &
    Matrix<double>::
    AddMMMul(const Matrix<double> & a, const Matrix<double>& b);

  template<>
    Matrix<float> &
    Matrix<float>::
    RepMMMul(const Matrix<float> & a, const Matrix<float>& b);

  template<>
    Matrix<float> &
    Matrix<float>::
    AddCMMtMul(const float& c, const Matrix<float>& a,
                               const Matrix<float>& b);

  template<>
    Matrix<double> &
    Matrix<double>::
    AddCMMtMul(const double& c, const Matrix<double>& a,
                                const Matrix<double>& b);

  template<>
    Matrix<float> &
    Matrix<float>::
    AddCVVtMul(const float& c, const BasicVector<float>& rA,
                               const BasicVector<float>& rB);

  template<>
    Matrix<double> &
    Matrix<double>::
    AddCVVtMul(const double& c, const BasicVector<double>& rA,
                                const BasicVector<double>& rB);

} // namespace STK

//#ifndef STK_Matrix_h
#endif

max_stars_repo_name: troylee/nnet-asr

max_issues_repo_name: troylee/nnet-asr

max_forks_repo_name: troylee/nnet-asr


text:
#include <mex.h>
#include <matrix.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <blas.h>    // dgemv
#include <lapack.h>  // dstegr

#define  abs1(a)         ((a) < 0.0 ? -(a) : (a))
#define  sign1(a)        ((a)==0) ? 0 : (((a)>0.0)?1:(-1))
#define  max1(a,b)       ((a) > (b) ? (a) : (b))
#define  min1(a,b)       ((a) < (b) ? (a) : (b))

double dot1(const double*a,const double*b,const ptrdiff_t n)
{
    double ret = 0;
    ptrdiff_t i;
    for(i=0;i<n;i++)
        ret +=a[i]*b[i];
    return ret;
}
void AxSym(const double*A, const double*x, double*result, const ptrdiff_t n)
{
    /*
     * input:
     * A: n x n
     * x: n x 1
     * output:
     * result = A*x: n x 1
     * NOTE: A is a symmetric matrix, only the the lower triangular part of A is to be referenced
     *       x and the result can not be the same space
     */
    char *UPLO="L";ptrdiff_t N;double ALPHA; ptrdiff_t LDA; ptrdiff_t INCX; double BETA;ptrdiff_t INCY;
    N = n; ALPHA = 1; LDA = n; INCX=1; BETA=0; INCY=1;
    dsymv(UPLO,&N,&ALPHA,A,&LDA,x,&INCX,&BETA,result,&INCY);
}

double computeObj(double s, double u0, double u1, double u2,double d0,double d1,double d2)
{
    return (u2*s*s + u1*s + u0)  /  (d2*s*s + d1*s + d0);
}
double quadfrac2(double u0, double u1, double u2,double d0,double d1,double d2)
{
    double x1;double x2;
    double c2 = u2*d1 - d2*u1;
    double c1 = 2*u2*d0 - 2*d2*u0;
    double c0 = u1*d0 - d1*u0;
    double dd = sqrt(c1*c1 - 4*c2*c0);

    if(c2==0)
    {
        x1 = -c0/c1;
        x2 = x1;
    }
    else
    {
        x1 = (-c1+dd)/(2*c2);
        x2 = (-c1-dd)/(2*c2);
    }
    if(computeObj(x1,u0,u1,u2,d0,d1,d2)<computeObj(x2,u0,u1,u2,d0,d1,d2))
    {
        return x1;
    }
    else
    {
        return x2;
    }

}

ptrdiff_t pos_max(const double*a,const ptrdiff_t n)
{
    double val = abs1(a[0]);
    ptrdiff_t pos=0;
    ptrdiff_t i ;
    for (i=1;i<n;i++)
    {
        if(abs1(a[i]) > val)
        {
            val = abs1(a[i]);
            pos = i;
        }
    }
    return pos;
}
void pvec(const double*a,const ptrdiff_t n)
{
    // print the vector a
    ptrdiff_t i;
    printf("[ ");
    for (i=0;i<n;i++)
    {
        printf("%f ",a[i]);
    }
    printf("]\n");
}

void solve(double *x,const double *in_x, const double*A, const double*b, const double c, const double *Q, const double *r,const double s,const ptrdiff_t n,const ptrdiff_t max_iter)
{

    double *Ax = (double*)malloc(n*sizeof(double));
    double *Qx = (double*)malloc(n*sizeof(double));
    double *up_g =  (double*)malloc(n*sizeof(double));
    double *down_g =  (double*)malloc(n*sizeof(double));
    double *grad =  (double*)malloc(n*sizeof(double));
    double xAx = 0;double bx = 0; double xQx = 0; double rx = 0;
	ptrdiff_t iter, i,j;

    memcpy(x,in_x,n*sizeof(double));
    AxSym(A, x, Ax, n);
    xAx = dot1(x,Ax,n);
    bx = dot1(x,b,n);
    AxSym(Q,x,Qx,n);
    xQx = dot1(x,Qx,n);
    rx = dot1(r,x,n);

    for ( iter=1;iter<=max_iter;iter++)
    {
        double u0,u1,u2,d0,d1,d2,alpha;
        double up,down,fobj;

        up = 0.5*xAx + bx + c;
        down = 0.5*xQx + rx + s;
        fobj =  up / down;

        for (j=0;j<n;j++)
        {
            up_g[j] = Ax[j]+b[j];
            down_g[j] = Qx[j]+r[j];
            grad[j] = up_g[j]/down - fobj * down_g[j]/down ;
        }
        i = pos_max(grad,n);
        u0 = 0.5*xAx + bx + c;
        u1 = Ax[i] + b[i];
        u2 = 0.5*A[i*n+i];
        d0 = 0.5*xQx + rx + s;
        d1 = Qx[i] + r[i];
        d2 = 0.5*Q[i*n+i];

        alpha =  quadfrac2(u0,u1,u2,d0,d1,d2);
        x[i] = x[i] + alpha;

        xAx = xAx + 2*alpha*Ax[i] + alpha*alpha*A[i*n+i];
        xQx = xQx + 2*alpha*Qx[i] + alpha*alpha*Q[i*n+i];
        for (j=0;j<n;j++)
        {
            Ax[j] = Ax[j] + alpha*A[i*n+j];
            Qx[j] = Qx[j] + alpha*Q[i*n+j];
        }
        bx = bx + alpha*b[i];
        rx = rx + alpha*r[i];

    }

    for (ptrdiff_t i=0;i<n;i++)if(x[i]==0)x[i]=1e-100;
    free(Ax);
    free(Qx);
    free(up_g);
    free(down_g);
    free(grad);

}

void mexFunction (int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
{
    /*set up input arguments */
    double* x =            mxGetPr(prhs[0]);
    double* A =            mxGetPr(prhs[1]);
    double* b =            mxGetPr(prhs[2]);
    double c =            mxGetScalar(prhs[3]);
    double* D =            mxGetPr(prhs[4]);
    double* e =            mxGetPr(prhs[5]);
    double f =            mxGetScalar(prhs[6]);
    ptrdiff_t n =         mxGetScalar(prhs[7]);
    ptrdiff_t max_iter =         mxGetScalar(prhs[8]);
    double *out_x;
    plhs[0] = mxCreateDoubleMatrix(n,1,mxREAL);
    out_x = mxGetPr(plhs[0]);
    solve(out_x,x,A,b,c,D,e,f,n,max_iter);
}

max_stars_repo_name: RVreugdenhil/sparseQP

max_issues_repo_name: RVreugdenhil/sparseQP

max_forks_repo_name: RVreugdenhil/sparseQP