GALAHAD SLLS package#
purpose#
The slls
package uses a preconditioned project-gradient method to solve
a given simplex-constrained linear least-squares problem.
The aim is to minimize the (regularized) least-squares objective function
See Section 4 of $GALAHAD/doc/slls.pdf for additional details.
terminology#
Any required solution
method#
Projected-gradient methods iterate towards a point that satisfies
these optimality conditions by ultimately aiming to satisfy
The method is iterative. Each iteration proceeds in two stages.
Firstly, a search direction
reference#
Full details are provided in
N. I. M. Gould (2023). ``Linear least-squares over the unit simplex’’. STFC-Rutherford Appleton Laboratory Computational Mathematics Group Internal Report 2023-2.
matrix storage#
The unsymmetric
Dense storage format:
The matrix
Dense by columns storage format:
The matrix
Sparse co-ordinate storage format:
Only the nonzero entries of the matrices are stored.
For the
Sparse row-wise storage format:
Again only the nonzero entries are stored, but this time
they are ordered so that those in row i appear directly before those
in row i+1. For the i-th row of
Sparse column-wise storage format:
Once again only the nonzero entries are stored, but this time
they are ordered so that those in column j appear directly before those
in column j+1. For the j-th column of
introduction to function calls#
To solve a given problem, functions from the slls package must be called in the following order:
slls_initialize - provide default control parameters and set up initial data structures
slls_read_specfile (optional) - override control values by reading replacement values from a file
set up problem data structures and fixed values by caling one of
slls_import - in the case that A
is explicitly availableslls_import_without_a - in the case that only the effect of applying A
and its transpose to a vector is possible
slls_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
solve the problem by calling one of
slls_solve_given_a - solve the problem using values of
slls_solve_reverse_a_prod - solve the problem by returning to the caller for products of A
and its transpose with specified vectors
slls_information (optional) - recover information about the solution and solution process
slls_terminate - deallocate data structures
See the examples section for illustrations of use.
callable functions#
overview of functions provided#
// typedefs typedef float spc_; typedef double rpc_; typedef int ipc_; // structs struct slls_control_type; struct slls_inform_type; struct slls_time_type; // function calls void slls_initialize( void **data, struct slls_control_type* control, ipc_ *status ); void slls_read_specfile( struct slls_control_type* control, const char specfile[] ); void slls_import( struct slls_control_type* control, void **data, ipc_ *status, ipc_ n, ipc_ o, const char A_type[], ipc_ Ao_ne, const ipc_ Ao_row[], const ipc_ Ao_col[], ipc_ Ao_ptr_ne, const ipc_ Ao_ptr[] ); void slls_import_without_a( struct slls_control_type* control, void **data, ipc_ *status, ipc_ n, ipc_ m ); void slls_reset_control( struct slls_control_type* control, void **data, ipc_ *status ); void slls_solve_given_a( void **data, void *userdata, ipc_ *status, ipc_ n, ipc_ o, ipc_ Ao_ne, const rpc_ Ao_val[], const rpc_ b[], rpc_ x[], rpc_ z[], rpc_ r[], rpc_ g[], ipc_ x_stat[], const rpc_ w[], ipc_(*)(ipc_, const rpc_[], rpc_[], const void*) eval_prec ); void slls_solve_reverse_a_prod( void **data, ipc_ *status, ipc_ *eval_status, ipc_ n, ipc_ o, const rpc_ b[], const rpc_ x_l[], const rpc_ x_u[], rpc_ x[], rpc_ z[], rpc_ r[], rpc_ g[], ipc_ x_stat[], rpc_ v[], const rpc_ p[], ipc_ nz_v[], ipc_ *nz_v_start, ipc_ *nz_v_end, const ipc_ nz_p[], ipc_ nz_p_end, const rpc_ w[] ); void slls_information(void **data, struct slls_inform_type* inform, ipc_ *status); void slls_terminate( void **data, struct slls_control_type* control, struct slls_inform_type* inform );
typedefs#
typedef float spc_
spc_
is real single precision
typedef double rpc_
rpc_
is the real working precision used, but may be changed to float
by
defining the preprocessor variable REAL_32
or (if supported) to
__real128
using the variable REAL_128
.
typedef int ipc_
ipc_
is the default integer word length used, but may be changed to
int64_t
by defining the preprocessor variable INTEGER_64
.
function calls#
void slls_initialize( void **data, struct slls_control_type* control, ipc_ *status )
Set default control values and initialize private data
Parameters:
data |
holds private internal data |
control |
is a struct containing control information (see slls_control_type) |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently):
|
void slls_read_specfile( struct slls_control_type* control, const char specfile[] )
Read the content of a specification file, and assign values associated with given keywords to the corresponding control parameters. An in-depth discussion of specification files is available, and a detailed list of keywords with associated default values is provided in $GALAHAD/src/slls/SLLS.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/slls.pdf for a list of how these keywords relate to the components of the control structure.
Parameters:
control |
is a struct containing control information (see slls_control_type) |
specfile |
is a character string containing the name of the specification file |
void slls_import( struct slls_control_type* control, void **data, ipc_ *status, ipc_ n, ipc_ o, const char Ao_type[], ipc_ Ao_ne, const ipc_ Ao_row[], const ipc_ Ao_col[], ipc_ Ao_ptr_ne, const ipc_ Ao_ptr[] )
Import problem data into internal storage prior to solution.
Parameters:
control |
is a struct whose members provide control paramters for the remaining prcedures (see slls_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
n |
is a scalar variable of type ipc_, that holds the number of variables. |
o |
is a scalar variable of type ipc_, that holds the number of residuals. |
Ao_type |
is a one-dimensional array of type char that specifies the symmetric storage scheme used for the design matrix |
Ao_ne |
is a scalar variable of type ipc_, that holds the number of entries in |
Ao_row |
is a one-dimensional array of size Ao_ne and type ipc_, that holds the row indices of |
Ao_col |
is a one-dimensional array of size Ao_ne and type ipc_, that holds the column indices of |
Ao_ptr_ne |
is a scalar variable of type ipc_, that holds the length of the pointer array if sparse row or column storage scheme is used for |
Ao_ptr |
is a one-dimensional array of size Ao_ptr_ne and type ipc_, that holds the starting position of each row of |
void slls_import_without_a( struct slls_control_type* control, void **data, ipc_ *status, ipc_ n, ipc_ o )
Import problem data into internal storage prior to solution.
Parameters:
control |
is a struct whose members provide control paramters for the remaining prcedures (see slls_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
n |
is a scalar variable of type ipc_, that holds the number of variables. |
o |
is a scalar variable of type ipc_, that holds the number of residuals. |
void slls_reset_control( struct slls_control_type* control, void **data, ipc_ *status )
Reset control parameters after import if required.
Parameters:
control |
is a struct whose members provide control paramters for the remaining prcedures (see slls_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
void slls_solve_given_a( void **data, void *userdata, ipc_ *status, ipc_ n, ipc_ o, ipc_ Ao_ne, const rpc_ Ao_val[], const rpc_ b[], rpc_ x[], rpc_ z[], rpc_ r[], rpc_ g[], ipc_ x_stat[], const rpc_ w[], ipc_(*)(ipc_, const rpc_[], rpc_[], const void*) eval_prec )
Solve the simplex-constrained linear least-squares problem when the design matrix
Parameters:
data |
holds private internal data |
userdata |
is a structure that allows data to be passed into the function and derivative evaluation programs. |
status |
is a scalar variable of type ipc_, that gives the entry and exit status from the package. On initial entry, status must be set to 1. Possible exit values are:
|
n |
is a scalar variable of type ipc_, that holds the number of variables |
o |
is a scalar variable of type ipc_, that holds the number of residuals. |
Ao_ne |
is a scalar variable of type ipc_, that holds the number of entries in the design matrix |
Ao_val |
is a one-dimensional array of size Ao_ne and type rpc_, that holds the values of the entries in the design matrix |
b |
is a one-dimensional array of size o and type rpc_, that holds the constant term |
x |
is a one-dimensional array of size n and type rpc_, that holds the values |
z |
is a one-dimensional array of size n and type rpc_, that holds the values |
r |
is a one-dimensional array of size o and type rpc_, that holds the values of the residuals |
g |
is a one-dimensional array of size n and type rpc_, that holds the values of the gradient |
x_stat |
is a one-dimensional array of size n and type ipc_, that gives the optimal status of the problem variables. If x_stat(j) is negative, the variable |
eval_prec |
is an optional user-supplied function that may be NULL. If non-NULL, it must have the following signature: ipc_ eval_prec( ipc_ n, const rpc_ v[], rpc_ p[], const void *userdata ) The product |
void slls_solve_reverse_a_prod( void **data, ipc_ *status, ipc_ *eval_status, ipc_ n, ipc_ o, const rpc_ b[], const rpc_ x_l[], const rpc_ x_u[], rpc_ x[], rpc_ z[], rpc_ r[], rpc_ g[], ipc_ x_stat[], rpc_ v[], const rpc_ p[], ipc_ nz_v[], ipc_ *nz_v_start, ipc_ *nz_v_end, const ipc_ nz_p[], ipc_ nz_p_end, const rpc_ w[] )
Solve the bound-constrained linear least-squares problem when the products of the Jacobian
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the entry and exit status from the package. Possible exit values are:
|
eval_status |
is a scalar variable of type ipc_, that is used to indicate if the matrix products can be provided (see |
n |
is a scalar variable of type ipc_, that holds the number of variables |
o |
is a scalar variable of type ipc_, that holds the number of residuals. |
b |
is a one-dimensional array of size o and type rpc_, that holds the constant term |
x_l |
is a one-dimensional array of size n and type rpc_, that holds the lower bounds |
x_u |
is a one-dimensional array of size n and type rpc_, that holds the upper bounds |
x |
is a one-dimensional array of size n and type rpc_, that holds the values |
r |
is a one-dimensional array of size m and type rpc_, that holds the values of the residuals |
g |
is a one-dimensional array of size n and type rpc_, that holds the values of the gradient |
z |
is a one-dimensional array of size n and type rpc_, that holds the values |
x_stat |
is a one-dimensional array of size n and type ipc_, that gives the optimal status of the problem variables. If x_stat(j) is negative, the variable |
v |
is a one-dimensional array of size n and type rpc_, that is used for reverse communication (see status=2-4 above for details). |
p |
is a one-dimensional array of size n and type rpc_, that is used for reverse communication (see status=2-4 above for details). |
nz_v |
is a one-dimensional array of size n and type ipc_, that is used for reverse communication (see status=3-4 above for details). |
nz_v_start |
is a scalar of type ipc_, that is used for reverse communication (see status=3-4 above for details). |
nz_v_end |
is a scalar of type ipc_, that is used for reverse communication (see status=3-4 above for details). |
nz_p |
is a one-dimensional array of size n and type ipc_, that is used for reverse communication (see status=4 above for details). |
nz_p_end |
is a scalar of type ipc_, that is used for reverse communication (see status=4 above for details). |
w |
is an optional one-dimensional array of size m and type rpc_, that holds the values |
void slls_information(void **data, struct slls_inform_type* inform, ipc_ *status)
Provides output information
Parameters:
data |
holds private internal data |
inform |
is a struct containing output information (see slls_inform_type) |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently):
|
void slls_terminate( void **data, struct slls_control_type* control, struct slls_inform_type* inform )
Deallocate all internal private storage
Parameters:
data |
holds private internal data |
control |
is a struct containing control information (see slls_control_type) |
inform |
is a struct containing output information (see slls_inform_type) |
available structures#
slls_control_type structure#
#include <galahad_slls.h> struct slls_control_type { // components bool f_indexing; ipc_ error; ipc_ out; ipc_ print_level; ipc_ start_print; ipc_ stop_print; ipc_ print_gap; ipc_ maxit; ipc_ cold_start; ipc_ preconditioner; ipc_ ratio_cg_vs_sd; ipc_ change_max; ipc_ cg_maxit; ipc_ arcsearch_max_steps; ipc_ sif_file_device; rpc_ weight; rpc_ stop_d; rpc_ stop_cg_relative; rpc_ stop_cg_absolute; rpc_ alpha_max; rpc_ alpha_initial; rpc_ alpha_reduction; rpc_ arcsearch_acceptance_tol; rpc_ stabilisation_weight; rpc_ cpu_time_limit; bool direct_subproblem_solve; bool exact_arc_search; bool space_critical; bool deallocate_error_fatal; bool generate_sif_file; char sif_file_name[31]; char prefix[31]; struct sbls_control_type sbls_control; struct convert_control_type convert_control; };
detailed documentation#
control derived type as a C struct
components#
bool f_indexing
use C or Fortran sparse matrix indexing
ipc_ error
unit number for error and warning diagnostics
ipc_ out
general output unit number
ipc_ print_level
the level of output required
ipc_ start_print
on which iteration to start printing
ipc_ stop_print
on which iteration to stop printing
ipc_ print_gap
how many iterations between printing
ipc_ maxit
how many iterations to perform (-ve reverts to HUGE(1)-1)
ipc_ cold_start
cold_start should be set to 0 if a warm start is required (with variable assigned according to X_stat, see below), and to any other value if the values given in prob.X suffice
ipc_ preconditioner
the preconditioner (scaling) used. Possible values are: /li 0. no preconditioner. /li 1. a diagonal preconditioner that normalizes the rows of
ipc_ ratio_cg_vs_sd
the ratio of how many iterations use CGLS rather than steepest descent
ipc_ change_max
the maximum number of per-iteration changes in the working set permitted when allowing CGLS rather than steepest descent
ipc_ cg_maxit
how many CG iterations to perform per SLLS iteration (-ve reverts to n+1)
ipc_ arcsearch_max_steps
the maximum number of steps allowed in a piecewise arcsearch (-ve=infini
ipc_ sif_file_device
the unit number to write generated SIF file describing the current probl
rpc_ weight
the value of the non-negative regularization weight sigma, i.e., the quadratic objective function q(x) will be regularized by adding 1/2 weight ||x||^2; any value smaller than zero will be regarded as zero.
rpc_ stop_d
the required accuracy for the dual infeasibility
rpc_ stop_cg_relative
the CG iteration will be stopped as soon as the current norm of the preconditioned gradient is smaller than max( stop_cg_relative * initial preconditioned gradient, stop_cg_absolute)
rpc_ alpha_max
the largest permitted arc length during the piecewise line search
rpc_ alpha_initial
the initial arc length during the inexact piecewise line search
rpc_ alpha_reduction
the arc length reduction factor for the inexact piecewise line search
rpc_ arcsearch_acceptance_tol
the required relative reduction during the inexact piecewise line search
rpc_ stabilisation_weight
the stabilisation weight added to the search-direction subproblem
rpc_ cpu_time_limit
the maximum CPU time allowed (-ve = no limit)
bool direct_subproblem_solve
direct_subproblem_solve is true if the least-squares subproblem is to be solved using a matrix factorization, and false if conjugate gradients are to be preferred
bool exact_arc_search
exact_arc_search is true if an exact arc_search is required, and false if an approximation suffices
bool space_critical
if space_critical is true, every effort will be made to use as little space as possible. This may result in longer computation times
bool deallocate_error_fatal
if deallocate_error_fatal is true, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue
bool generate_sif_file
if generate_sif_file is true, a SIF file describing the current problem will be generated
char sif_file_name[31]
name (max 30 characters) of generated SIF file containing input problem
char prefix[31]
all output lines will be prefixed by a string (max 30 characters) prefix(2:LEN(TRIM(.prefix))-1) where prefix contains the required string enclosed in quotes, e.g. “string” or ‘string’
struct sbls_control_type sbls_control
control parameters for SBLS
struct convert_control_type convert_control
control parameters for CONVERT
slls_time_type structure#
#include <galahad_slls.h> struct slls_time_type { // components rpc_ total; rpc_ analyse; rpc_ factorize; rpc_ solve; };
detailed documentation#
time derived type as a C struct
components#
rpc_ total
the total CPU time spent in the package
rpc_ analyse
the CPU time spent analysing the required matrices prior to factorization
rpc_ factorize
the CPU time spent factorizing the required matrices
rpc_ solve
the CPU time spent in the linear solution phase
slls_inform_type structure#
#include <galahad_slls.h> struct slls_inform_type { // components ipc_ status; ipc_ alloc_status; ipc_ factorization_status; ipc_ iter; ipc_ cg_iter; rpc_ obj; rpc_ norm_pg; char bad_alloc[81]; struct slls_time_type time; struct sbls_inform_type sbls_inform; struct convert_inform_type convert_inform; };
detailed documentation#
inform derived type as a C struct
components#
ipc_ status
reported return status.
ipc_ alloc_status
Fortran STAT value after allocate failure.
ipc_ factorization_status
status return from factorization
ipc_ iter
number of iterations required
ipc_ cg_iter
number of CG iterations required
rpc_ obj
current value of the objective function, r(x).
rpc_ norm_pg
current value of the Euclidean norm of projected gradient of r(x).
char bad_alloc[81]
name of array which provoked an allocate failure
struct slls_time_type time
times for various stages
struct sbls_inform_type sbls_inform
inform values from SBLS
struct convert_inform_type convert_inform
inform values for CONVERT
example calls#
This is an example of how to use the package to solve a bound-constrained linear least-squares problem; the code is available in $GALAHAD/src/slls/C/sllst.c . A variety of supported Hessian and constraint matrix storage formats are shown.
Notice that C-style indexing is used, and that this is flagged by setting
control.f_indexing
to false
. The floating-point type rpc_
is set in galahad_precision.h
to double
by default, but to float
if the preprocessor variable SINGLE
is defined. Similarly, the integer
type ipc_
from galahad_precision.h
is set to int
by default,
but to int64_t
if the preprocessor variable INTEGER_64
is defined.
/* sllst.c */
/* Full test for the SLLS C interface using C sparse matrix indexing */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_slls.h"
#ifdef REAL_128
#include <quadmath.h>
#endif
// Define imax
ipc_ imax(ipc_ a, ipc_ b) {
return (a > b) ? a : b;
};
// Custom userdata struct
struct userdata_type {
rpc_ scale;
};
// Function prototypes
ipc_ prec( ipc_ n, const rpc_ v[], rpc_ p[], const void * );
int main(void) {
// Derived types
void *data;
struct slls_control_type control;
struct slls_inform_type inform;
// Set user data
struct userdata_type userdata;
userdata.scale = 1.0;
// Set problem data
ipc_ n = 10; // dimension
ipc_ o = n + 1; // number of residuals
ipc_ Ao_ne = 2 * n; // sparse Jacobian elements
ipc_ Ao_dense_ne = o * n; // dense Jacobian elements
// row-wise storage
ipc_ Ao_row[Ao_ne]; // row indices,
ipc_ Ao_col[Ao_ne]; // column indices
ipc_ Ao_ptr_ne = o+1; // row pointer length
ipc_ Ao_ptr[Ao_ptr_ne]; // row pointers
rpc_ Ao_val[Ao_ne]; // values
rpc_ Ao_dense[Ao_dense_ne]; // dense values
// column-wise storage
ipc_ Ao_by_col_row[Ao_ne]; // row indices,
ipc_ Ao_by_col_ptr_ne = n+1; // column pointer length
ipc_ Ao_by_col_ptr[Ao_by_col_ptr_ne]; // column pointers
rpc_ Ao_by_col_val[Ao_ne]; // values
rpc_ Ao_by_col_dense[Ao_dense_ne]; // dense values
rpc_ b[o]; // linear term in the objective
rpc_ x[n]; // variables
rpc_ z[n]; // dual variables
rpc_ r[o]; // residual
rpc_ g[n]; // gradient
// Set output storage
ipc_ x_stat[n]; // variable status
char st[3];
ipc_ i, l, status;
// A = ( I ) and b = ( i * e )
// ( e^T ) ( n + 1 )
for( ipc_ i = 0; i < n; i++) b[i] = i + 1;
b[n] = n+1;
// A by rows
for( ipc_ i = 0; i < n; i++)
{
Ao_ptr[i] = i;
Ao_row[i] = i; Ao_col[i] = i; Ao_val[i] = 1.0;
}
Ao_ptr[n] = n;
for( ipc_ i = 0; i < n; i++)
{
Ao_row[n+i] = n; Ao_col[n+i] = i; Ao_val[n+i] = 1.0;
}
Ao_ptr[o] = Ao_ne;
l = - 1;
for( ipc_ i = 0; i < n; i++)
{
for( ipc_ j = 0; j < n; j++)
{
l = l + 1;
if ( i == j ) {
Ao_dense[l] = 1.0;
}
else {
Ao_dense[l] = 0.0;
}
}
}
for( ipc_ j = 0; j < n; j++)
{
l = l + 1;
Ao_dense[l] = 1.0;
}
// A by columns
l = - 1;
for( ipc_ j = 0; j < n; j++)
{
l = l + 1; Ao_by_col_ptr[j] = l ;
Ao_by_col_row[l] = j ; Ao_by_col_val[l] = 1.0;
l = l + 1;
Ao_by_col_row[l] = n ; Ao_by_col_val[l] = 1.0;
}
Ao_by_col_ptr[n] = Ao_ne;
l = - 1;
for( ipc_ j = 0; j < n; j++)
{
for( ipc_ i = 0; i < n; i++)
{
l = l + 1;
if ( i == j ) {
Ao_by_col_dense[l] = 1.0;
}
else {
Ao_by_col_dense[l] = 0.0;
}
}
l = l + 1;
Ao_by_col_dense[l] = 1.0;
}
printf(" C sparse matrix indexing\n\n");
printf(" basic tests of slls storage formats\n\n");
for( ipc_ d=1; d <= 5; d++){
// Initialize SLLS
slls_initialize( &data, &control, &status );
// Set user-defined control options
control.f_indexing = false; // C sparse matrix indexing
// Start from 0
for( ipc_ i = 0; i < n; i++) x[i] = 0.0;
for( ipc_ i = 0; i < n; i++) z[i] = 0.0;
switch(d){
case 1: // sparse co-ordinate storage
strcpy( st, "CO" );
slls_import( &control, &data, &status, n, o,
"coordinate", Ao_ne, Ao_row, Ao_col, 0, NULL );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_ne, Ao_val, b,
x, z, r, g, x_stat, prec );
break;
case 2: // sparse by rows
strcpy( st, "SR" );
slls_import( &control, &data, &status, n, o,
"sparse_by_rows", Ao_ne, NULL, Ao_col,
Ao_ptr_ne, Ao_ptr );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_ne, Ao_val, b,
x, z, r, g, x_stat, prec );
break;
case 3: // dense by rows
strcpy( st, "DR" );
slls_import( &control, &data, &status, n, o,
"dense_by_rows", Ao_dense_ne,
NULL, NULL, 0, NULL );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_dense_ne, Ao_dense, b,
x, z, r, g, x_stat, prec );
break;
case 4: // sparse by columns
strcpy( st, "SC" );
slls_import( &control, &data, &status, n, o,
"sparse_by_columns", Ao_ne, Ao_by_col_row,
NULL, Ao_by_col_ptr_ne, Ao_by_col_ptr );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_ne, Ao_by_col_val, b,
x, z, r, g, x_stat, prec );
break;
case 5: // dense by columns
strcpy( st, "DC" );
slls_import( &control, &data, &status, n, o,
"dense_by_columns", Ao_dense_ne,
NULL, NULL, 0, NULL );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_dense_ne, Ao_by_col_dense, b,
x, z, r, g, x_stat, prec );
break;
}
slls_information( &data, &inform, &status );
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_sf.h
#include "galahad_pquad_sf.h"
#else
printf("%s:%6" i_ipc_ " iterations. Optimal objective value = %5.2f"
" status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%s: SLLS_solve exit status = %1" i_ipc_ "\n",
st, inform.status);
}
//printf("x: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);
//printf("\n");
//printf("gradient: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", g[i]);
//printf("\n");
// Delete internal workspace
slls_terminate( &data, &control, &inform );
}
printf("\n tests reverse-communication options\n\n");
// reverse-communication input/output
ipc_ on;
on = imax( n, o );
ipc_ eval_status, nz_v_start, nz_v_end, nz_p_end;
ipc_ nz_v[on], nz_p[o], mask[o];
rpc_ v[on], p[on];
nz_p_end = 0;
// Initialize SLLS
slls_initialize( &data, &control, &status );
// Set user-defined control options
control.f_indexing = false; // C sparse matrix indexing
// Start from 0
for( ipc_ i = 0; i < n; i++) x[i] = 0.0;
for( ipc_ i = 0; i < n; i++) z[i] = 0.0;
strcpy( st, "RC" );
for( ipc_ i = 0; i < o; i++) mask[i] = 0;
slls_import_without_a( &control, &data, &status, n, o ) ;
while(true){ // reverse-communication loop
slls_solve_reverse_a_prod( &data, &status, &eval_status, n, o, b,
x, z, r, g, x_stat, v, p,
nz_v, &nz_v_start, &nz_v_end,
nz_p, nz_p_end );
if(status == 0){ // successful termination
break;
}else if(status < 0){ // error exit
break;
}else if(status == 2){ // evaluate p = Av
p[n]=0.0;
for( ipc_ i = 0; i < n; i++){
p[i] = v[i];
p[n] = p[n] + v[i];
}
}else if(status == 3){ // evaluate p = A^Tv
for( ipc_ i = 0; i < n; i++) p[i] = v[i] + v[n];
}else if(status == 4){ // evaluate p = Av for sparse v
p[n]=0.0;
for( ipc_ i = 0; i < n; i++) p[i] = 0.0;
for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
i = nz_v[l];
p[i] = v[i];
p[n] = p[n] + v[i];
}
}else if(status == 5){ // evaluate p = sparse Av for sparse v
nz_p_end = 0;
for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
i = nz_v[l];
if (mask[i] == 0){
mask[i] = 1;
nz_p[nz_p_end] = i;
nz_p_end = nz_p_end + 1;
p[i] = v[i];
}
if (mask[n] == 0){
mask[n] = 1;
nz_p[nz_p_end] = n;
nz_p_end = nz_p_end + 1;
p[n] = v[i];
}else{
p[n] = p[n] + v[i];
}
}
for( ipc_ l = 0; l < nz_p_end; l++) mask[nz_p[l]] = 0;
}else if(status == 6){ // evaluate p = sparse A^Tv
for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
i = nz_v[l];
p[i] = v[i] + v[n];
}
}else if(status == 7){ // evaluate p = P^{-}v
for( ipc_ i = 0; i < n; i++) p[i] = userdata.scale * v[i];
}else{
printf(" the value %1" i_ipc_ " of status should not occur\n", status);
break;
}
eval_status = 0;
}
// Record solution information
slls_information( &data, &inform, &status );
// Print solution details
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_sf.h
#include "galahad_pquad_sf.h"
#else
printf("%s:%6" i_ipc_ " iterations. Optimal objective value = %5.2f"
" status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%s: SLLS_solve exit status = %1" i_ipc_ "\n",
st, inform.status);
}
//printf("x: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);
//printf("\n");
//printf("gradient: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", g[i]);
//printf("\n");
// Delete internal workspace
slls_terminate( &data, &control, &inform );
}
// Apply preconditioner
ipc_ prec( ipc_ n, const rpc_ v[], rpc_ p[], const void *userdata ){
struct userdata_type *myuserdata = (struct userdata_type *) userdata;
rpc_ scale = myuserdata->scale;
for( ipc_ i = 0; i < n; i++) p[i] = scale * v[i];
return 0;
}
This is the same example, but now fortran-style indexing is used; the code is available in $GALAHAD/src/slls/C/sllstf.c .
/* sllstf.c */
/* Full test for the SLLS C interface using fortran sparse matrix indexing */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_slls.h"
#ifdef REAL_128
#include <quadmath.h>
#endif
// Define imax
ipc_ imax(ipc_ a, ipc_ b) {
return (a > b) ? a : b;
};
// Custom userdata struct
struct userdata_type {
rpc_ scale;
};
// Function prototypes
ipc_ prec( ipc_ n, const rpc_ v[], rpc_ p[], const void * );
int main(void) {
// Derived types
void *data;
struct slls_control_type control;
struct slls_inform_type inform;
// Set user data
struct userdata_type userdata;
userdata.scale = 1.0;
// Set problem data
ipc_ n = 10; // dimension
ipc_ o = n + 1; // number of residuals
ipc_ Ao_ne = 2 * n; // sparse Jacobian elements
ipc_ Ao_dense_ne = o * n; // dense Jacobian elements
// row-wise storage
ipc_ Ao_row[Ao_ne]; // row indices,
ipc_ Ao_col[Ao_ne]; // column indices
ipc_ Ao_ptr_ne = o+1; // row pointer length
ipc_ Ao_ptr[Ao_ptr_ne]; // row pointers
rpc_ Ao_val[Ao_ne]; // values
rpc_ Ao_dense[Ao_dense_ne]; // dense values
// column-wise storage
ipc_ Ao_by_col_row[Ao_ne]; // row indices,
ipc_ Ao_by_col_ptr_ne = n+1; // column pointer length
ipc_ Ao_by_col_ptr[Ao_by_col_ptr_ne]; // column pointers
rpc_ Ao_by_col_val[Ao_ne]; // values
rpc_ Ao_by_col_dense[Ao_dense_ne]; // dense values
rpc_ b[o]; // linear term in the objective
rpc_ x[n]; // variables
rpc_ z[n]; // dual variables
rpc_ r[o]; // residual
rpc_ g[n]; // gradient
// Set output storage
ipc_ x_stat[n]; // variable status
char st[3];
ipc_ i, l, status;
// A = ( I ) and b = ( i * e )
// ( e^T ) ( n + 1 )
for( ipc_ i = 0; i < n; i++) b[i] = i + 1;
b[n] = n+1;
// A by rows
for( ipc_ i = 0; i < n; i++)
{
Ao_ptr[i] = i + 1;
Ao_row[i] = i + 1; Ao_col[i] = i + 1; Ao_val[i] = 1.0;
}
Ao_ptr[n] = n + 1;
for( ipc_ i = 0; i < n; i++)
{
Ao_row[n+i] = o; Ao_col[n+i] = i + 1; Ao_val[n+i] = 1.0;
}
Ao_ptr[o] = Ao_ne + 1;
l = - 1;
for( ipc_ i = 0; i < n; i++)
{
for( ipc_ j = 0; j < n; j++)
{
l = l + 1;
if ( i == j ) {
Ao_dense[l] = 1.0;
}
else {
Ao_dense[l] = 0.0;
}
}
}
for( ipc_ j = 0; j < n; j++)
{
l = l + 1;
Ao_dense[l] = 1.0;
}
// A by columns
l = - 1;
for( ipc_ j = 0; j < n; j++)
{
l = l + 1; Ao_by_col_ptr[j] = l + 1;
Ao_by_col_row[l] = j + 1; Ao_by_col_val[l] = 1.0;
l = l + 1;
Ao_by_col_row[l] = o; Ao_by_col_val[l] = 1.0;
}
Ao_by_col_ptr[n] = Ao_ne;
l = - 1;
for( ipc_ j = 0; j < n; j++)
{
for( ipc_ i = 0; i < n; i++)
{
l = l + 1;
if ( i == j ) {
Ao_by_col_dense[l] = 1.0;
}
else {
Ao_by_col_dense[l] = 0.0;
}
}
l = l + 1;
Ao_by_col_dense[l] = 1.0;
}
printf(" fortran sparse matrix indexing\n\n");
printf(" basic tests of slls storage formats\n\n");
for( ipc_ d=1; d <= 5; d++){
// Initialize SLLS
slls_initialize( &data, &control, &status );
// Set user-defined control options
control.f_indexing = true; // fortran sparse matrix indexing
// Start from 0
for( ipc_ i = 0; i < n; i++) x[i] = 0.0;
for( ipc_ i = 0; i < n; i++) z[i] = 0.0;
switch(d){
case 1: // sparse co-ordinate storage
strcpy( st, "CO" );
slls_import( &control, &data, &status, n, o,
"coordinate", Ao_ne, Ao_row, Ao_col, 0, NULL );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_ne, Ao_val, b,
x, z, r, g, x_stat, prec );
break;
case 2: // sparse by rows
strcpy( st, "SR" );
slls_import( &control, &data, &status, n, o,
"sparse_by_rows", Ao_ne, NULL, Ao_col,
Ao_ptr_ne, Ao_ptr );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_ne, Ao_val, b,
x, z, r, g, x_stat, prec );
break;
case 3: // dense by rows
strcpy( st, "DR" );
slls_import( &control, &data, &status, n, o,
"dense_by_rows", Ao_dense_ne,
NULL, NULL, 0, NULL );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_dense_ne, Ao_dense, b,
x, z, r, g, x_stat, prec );
break;
case 4: // sparse by columns
strcpy( st, "SC" );
slls_import( &control, &data, &status, n, o,
"sparse_by_columns", Ao_ne, Ao_by_col_row,
NULL, Ao_by_col_ptr_ne, Ao_by_col_ptr );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_ne, Ao_by_col_val, b,
x, z, r, g, x_stat, prec );
break;
case 5: // dense by columns
strcpy( st, "DC" );
slls_import( &control, &data, &status, n, o,
"dense_by_columns", Ao_dense_ne,
NULL, NULL, 0, NULL );
slls_solve_given_a( &data, &userdata, &status, n, o,
Ao_dense_ne, Ao_by_col_dense, b,
x, z, r, g, x_stat, prec );
break;
}
slls_information( &data, &inform, &status );
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_sf.h
#include "galahad_pquad_sf.h"
#else
printf("%s:%6" i_ipc_ " iterations. Optimal objective value = %.2f"
" status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%s: SLLS_solve exit status = %1" i_ipc_ "\n",
st, inform.status);
}
//printf("x: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);
//printf("\n");
//printf("gradient: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", g[i]);
//printf("\n");
// Delete internal workspace
slls_terminate( &data, &control, &inform );
}
printf("\n tests reverse-communication options\n\n");
// reverse-communication input/output
ipc_ on;
on = imax( o, n );
ipc_ eval_status, nz_v_start, nz_v_end, nz_p_end;
ipc_ nz_v[on], nz_p[o], mask[o];
rpc_ v[on], p[on];
nz_p_end = 0;
// Initialize SLLS
slls_initialize( &data, &control, &status );
// Set user-defined control options
control.f_indexing = true; // fortran sparse matrix indexing
// Start from 0
for( ipc_ i = 0; i < n; i++) x[i] = 0.0;
for( ipc_ i = 0; i < n; i++) z[i] = 0.0;
strcpy( st, "RC" );
for( ipc_ i = 0; i < o; i++) mask[i] = 0;
slls_import_without_a( &control, &data, &status, n, o ) ;
while(true){ // reverse-communication loop
slls_solve_reverse_a_prod( &data, &status, &eval_status, n, o, b,
x, z, r, g, x_stat, v, p,
nz_v, &nz_v_start, &nz_v_end,
nz_p, nz_p_end );
if(status == 0){ // successful termination
break;
}else if(status < 0){ // error exit
break;
}else if(status == 2){ // evaluate p = Av
p[n]=0.0;
for( ipc_ i = 0; i < n; i++){
p[i] = v[i];
p[n] = p[n] + v[i];
}
}else if(status == 3){ // evaluate p = A^Tv
for( ipc_ i = 0; i < n; i++) p[i] = v[i] + v[n];
}else if(status == 4){ // evaluate p = Av for sparse v
p[n]=0.0;
for( ipc_ i = 0; i < n; i++) p[i] = 0.0;
for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
i = nz_v[l]-1;
p[i] = v[i];
p[n] = p[n] + v[i];
}
}else if(status == 5){ // evaluate p = sparse Av for sparse v
nz_p_end = 0;
for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
i = nz_v[l]-1;
if (mask[i] == 0){
mask[i] = 1;
nz_p[nz_p_end] = i+1;
nz_p_end = nz_p_end + 1;
p[i] = v[i];
}
if (mask[n] == 0){
mask[n] = 1;
nz_p[nz_p_end] = o;
nz_p_end = nz_p_end + 1;
p[n] = v[i];
}else{
p[n] = p[n] + v[i];
}
}
for( ipc_ l = 0; l < nz_p_end; l++) mask[nz_p[l]] = 0;
}else if(status == 6){ // evaluate p = sparse A^Tv
for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
i = nz_v[l]-1;
p[i] = v[i] + v[n];
}
}else if(status == 7){ // evaluate p = P^{-}v
for( ipc_ i = 0; i < n; i++) p[i] = userdata.scale * v[i];
}else{
printf(" the value %1" i_ipc_ " of status should not occur\n", status);
break;
}
eval_status = 0;
}
// Record solution information
slls_information( &data, &inform, &status );
// Print solution details
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_sf.h
#include "galahad_pquad_sf.h"
#else
printf("%s:%6" i_ipc_ " iterations. Optimal objective value = %.2f"
" status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%s: SLLS_solve exit status = %1" i_ipc_ "\n",
st, inform.status);
}
//printf("x: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);
//printf("\n");
//printf("gradient: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", g[i]);
//printf("\n");
// Delete internal workspace
slls_terminate( &data, &control, &inform );
}
// Apply preconditioner
ipc_ prec( ipc_ n, const rpc_ v[], rpc_ p[], const void *userdata ){
struct userdata_type *myuserdata = (struct userdata_type *) userdata;
rpc_ scale = myuserdata->scale;
for( ipc_ i = 0; i < n; i++) p[i] = scale * v[i];
return 0;
}