GALAHAD

https://ralna.github.io/galahad_docs/
Admin
GALAHAD is an open-source library of thread-safe packages for solving continuous nonlinear...
Tokens:2,227,369
Snippets:14,161
Trust Score:9.7
Update:1 month ago
Show doc for...
Tokens:
Raw
### DGO C interface example initialization

Source: https://ralna.github.io/galahad_docs/html/C/dgo.html

Example setup for the DGO C interface, including custom userdata structure definition.

```c
/* dgot.c */
/* Full test for the DGO C interface using C sparse matrix indexing */
/* Jari Fowkes & Nick Gould, STFC-Rutherford Appleton Laboratory, 2021 */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_dgo.h"
#ifdef REAL_128
#include <quadmath.h>
#endif

// Custom userdata struct
struct userdata_type {
   rpc_ p;
   rpc_ freq;
   rpc_ mag;
};
```

--------------------------------

### DGO C interface setup with Fortran indexing

Source: https://ralna.github.io/galahad_docs/html/C/dgo.html

Example setup for the DGO C interface using Fortran-style sparse matrix indexing, located in dgotf.c.

```C
/* dgotf.c */
/* Full test for the DGO C interface using Fortran sparse matrix indexing */
/* Jari Fowkes & Nick Gould, STFC-Rutherford Appleton Laboratory, 2021 */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_dgo.h"
#ifdef REAL_128
#include <quadmath.h>
#endif

// Custom userdata struct
struct userdata_type {
   rpc_ p;
   rpc_ freq;
   rpc_ mag;
};
```

--------------------------------

### DQP example calls: switch over storage schemes and shifted least-distance test setup

Source: https://ralna.github.io/galahad_docs/html/C/dqp.html

The first block shows the main DQP example loop that selects among storage schemes (coordinate, sparse_by_rows, dense, diagonal, scaled_identity, identity) via a switch on d, calls dqp_import and dqp_solve_qp for each, then reports results with dqp_information and cleans up with dqp_terminate. The second block starts the separate shifted least-distance test: it initializes DQP, sets control options (including f_indexing=false and sytr solvers), and initializes x, y, z arrays to zero before setting shifted least-distance data.

```c
switch(d){
            case 1: // sparse co-ordinate storage
                st = 'C';
                dqp_import( &control, &data, &status, n, m,
                           "coordinate", H_ne, H_row, H_col, NULL,
                           "coordinate", A_ne, A_row, A_col, NULL );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;
            printf(" case %1" d_ipc_ " break\n",d);
            case 2: // sparse by rows
                st = 'R';
                dqp_import( &control, &data, &status, n, m,
                            "sparse_by_rows", H_ne, NULL, H_col, H_ptr,
                            "sparse_by_rows", A_ne, NULL, A_col, A_ptr );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;
            case 3: // dense
                st = 'D';
                ipc_ H_dense_ne = 6; // number of elements of H
                ipc_ A_dense_ne = 6; // number of elements of A
                rpc_ H_dense[] = {1.0, 0.0, 1.0, 0.0, 0.0, 1.0};
                rpc_ A_dense[] = {2.0, 1.0, 0.0, 0.0, 1.0, 1.0};
                dqp_import( &control, &data, &status, n, m,
                            "dense", H_ne, NULL, NULL, NULL,
                            "dense", A_ne, NULL, NULL, NULL );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u,
                              x, c, y, z, x_stat, c_stat );
                break;
            case 4: // diagonal
                st = 'L';
                dqp_import( &control, &data, &status, n, m,
                            "diagonal", H_ne, NULL, NULL, NULL,
                            "sparse_by_rows", A_ne, NULL, A_col, A_ptr );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;

            case 5: // scaled identity
                st = 'S';
                dqp_import( &control, &data, &status, n, m,
                            "scaled_identity", H_ne, NULL, NULL, NULL,
                            "sparse_by_rows", A_ne, NULL, A_col, A_ptr );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;
            case 6: // identity
                st = 'I';
                dqp_import( &control, &data, &status, n, m,
                            "identity", H_ne, NULL, NULL, NULL,
                            "sparse_by_rows", A_ne, NULL, A_col, A_ptr );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;
            }
        dqp_information( &data, &inform, &status );

        if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_f.h
#include "galahad_pquad_f.h"
#else
            printf("%c:%6" d_ipc_ " iterations. Optimal objective "
                   "value = %.2f status = %1" d_ipc_ "\n",
                   st, inform.iter, inform.obj, inform.status);
#endif
        }else{
            printf("%c: DQP_solve exit status = %1" d_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
        dqp_terminate( &data, &control, &inform );
    }
```

```c
printf("\n separate test of sldqp\n\n");

    // test shifted least-distance interface
//  for( ipc_ d=1; d <= 0; d++){
    for( ipc_ d=1; d <= 1; d++){

        // Initialize DQP
        dqp_initialize( &data, &control, &status );

        // Set user-defined control options
        control.f_indexing = false; // C sparse matrix indexing
        strcpy(control.symmetric_linear_solver, "sytr ") ;
        strcpy(control.definite_linear_solver, "sytr ") ;
        strcpy(control.fdc_control.symmetric_linear_solver, "sytr ") ;
        control.fdc_control.use_sls = true;
        strcpy(control.sbls_control.symmetric_linear_solver, "sytr ") ;
        strcpy(control.sbls_control.definite_linear_solver, "sytr ") ;

        // Start from 0
        rpc_ x[] = {0.0,0.0,0.0};
        rpc_ y[] = {0.0,0.0};
        rpc_ z[] = {0.0,0.0,0.0};

        // Set shifted least-distance data
```

--------------------------------

### BQP Bound-Constrained QP Solver Example

Source: https://ralna.github.io/galahad_docs/html/Julia/bqp.html

A complete example demonstrating the setup and solution of a bound-constrained quadratic programming problem using different matrix storage formats and reverse-communication.

```julia
# test_bqp.jl
# Simple code to test the Julia interface to BQP

using GALAHAD
using Test
using Printf
using Accessors
using Quadmath

function test_bqp(::Type{T}, ::Type{INT}) where {T,INT}
  # Derived types
  data = Ref{Ptr{Cvoid}}()
  control = Ref{bqp_control_type{T,INT}}()
  inform = Ref{bqp_inform_type{T,INT}}()

  # Set problem data
  n = INT(10)  # dimension
  H_ne = 2 * n - 1 # Hesssian elements, NB lower triangle
  H_dense_ne = div(n * (n + 1), 2) # dense Hessian elements
  H_row = zeros(INT, H_ne) # row indices,
  H_col = zeros(INT, H_ne) # column indices
  H_ptr = zeros(INT, n + 1)  # row pointers
  H_val = zeros(T, H_ne) # values
  H_dense = zeros(T, H_dense_ne) # dense values
  H_diag = zeros(T, n)   # diagonal values
  g = zeros(T, n)  # linear term in the objective
  f = one(T)  # constant term in the objective
  x_l = zeros(T, n) # variable lower bound
  x_u = zeros(T, n) # variable upper bound
  x = zeros(T, n) # variables
  z = zeros(T, n) # dual variables

  # Set output storage
  x_stat = zeros(INT, n) # variable status
  st = ' '
  status = Ref{INT}()

  g[1] = 2.0
  for i in 2:n
    g[i] = 0.0
  end
  x_l[1] = -1.0
  for i in 2:n
    x_l[i] = -Inf
  end
  x_u[1] = 1.0
  x_u[2] = Inf
  for i in 3:n
    x_u[i] = 2.0
  end

  # H = tridiag(2,1), H_dense = diag(2)
  l = 1
  H_ptr[1] = l
  H_row[l] = 1
  H_col[l] = 1
  H_val[l] = 2.0
  for i in 2:n
    l = l + 1
    H_ptr[i] = l
    H_row[l] = i
    H_col[l] = i - 1
    H_val[l] = 1.0
    l = l + 1
    H_row[l] = i
    H_col[l] = i
    H_val[l] = 2.0
  end
  H_ptr[n + 1] = l + 1

  l = 0
  for i in 1:n
    H_diag[i] = 2.0
    for j in 1:i
      l = l + 1
      if j < i - 1
        H_dense[l] = 0.0
      elseif j == i - 1
        H_dense[l] = 1.0
      else
        H_dense[l] = 2.0
      end
    end
  end

  @printf(" fortran sparse matrix indexing\n\n")
  @printf(" basic tests of bqp storage formats\n\n")

  for d in 1:4
    # Initialize BQP
    bqp_initialize(T, INT, data, control, status)

    # Start from 0
    for i in 1:n
      x[i] = 0.0
      z[i] = 0.0
    end

    # sparse co-ordinate storage
    if d == 1
      st = 'C'
      bqp_import(T, INT, control, data, status, n,
                 "coordinate", H_ne, H_row, H_col, C_NULL)

      bqp_solve_given_h(T, INT, data, status, n, H_ne, H_val, g, f,
                        x_l, x_u, x, z, x_stat)
    end

    # sparse by rows
    if d == 2
      st = 'R'
      bqp_import(T, INT, control, data, status, n,
                 "sparse_by_rows", H_ne, C_NULL, H_col, H_ptr)

      bqp_solve_given_h(T, INT, data, status, n, H_ne, H_val, g, f,
                        x_l, x_u, x, z, x_stat)
    end

    # dense
    if d == 3
      st = 'D'
      bqp_import(T, INT, control, data, status, n,
                 "dense", H_dense_ne, C_NULL, C_NULL, C_NULL)

      bqp_solve_given_h(T, INT, data, status, n, H_dense_ne, H_dense,
                        g, f, x_l, x_u, x, z, x_stat)
    end

    # diagonal
    if d == 4
      st = 'L'
      bqp_import(T, INT, control, data, status, n,
                 "diagonal", H_ne, C_NULL, C_NULL, C_NULL)

      bqp_solve_given_h(T, INT, data, status, n, n, H_diag, g, f,
                        x_l, x_u, x, z, x_stat)
    end

    bqp_information(T, INT, data, inform, status)

    if inform[].status == 0
      @printf("%c:%6i iterations. Optimal objective value = %5.2f status = %1i\n", st,
              inform[].iter, inform[].obj, inform[].status)
    else
      @printf("%c: BQP_solve exit status = %1i\n", st, inform[].status)
    end

    # Delete internal workspace
    bqp_terminate(T, INT, data, control, inform)
  end

  @printf("\n tests reverse-communication options\n\n")

  # reverse-communication input/output
  nz_v_start = Ref{INT}()
  nz_v_end = Ref{INT}()
  nz_v = zeros(INT, n)
  nz_prod = zeros(INT, n)
  mask = zeros(INT, n)
  v = zeros(T, n)
  prod = zeros(T, n)

  nz_prod_end = 1

  # Initialize BQP
  bqp_initialize(T, INT, data, control, status)
  # @reset control[].print_level = 1

  # Start from 0
  for i in 1:n
    x[i] = 0.0
    z[i] = 0.0
  end

  st = 'I'
  for i in 1:n
    mask[i] = 0
  end
  bqp_import_without_h(T, INT, control, data, status, n)

  terminated = false
  while !terminated # reverse-communication loop
    bqp_solve_reverse_h_prod(T, INT, data, status, n, g, f, x_l, x_u, x, z, x_stat, v, prod, nz_v,
                             nz_v_start, nz_v_end, nz_prod, nz_prod_end)
```

--------------------------------

### SBLS C interface example: sblst.c setup

Source: https://ralna.github.io/galahad_docs/html/C/sbls.html

This is the beginning of the example program sblst.c that solves a block system of linear equations using the SBLS C interface. It sets up problem data in various sparse and dense formats, initializes SBLS, and sets control options including C-style indexing (control.f_indexing = false) and the use of the sytr solvers.

```c
/* sblst.c */
/* Full test for the SBLS 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_sbls.h"
#ifdef REAL_128
#include <quadmath.h>
#endif

int main(void) {

    // Derived types
    void *data;
    struct sbls_control_type control;
    struct sbls_inform_type inform;

    // Set problem data
    ipc_ n = 3; // dimension of H
    ipc_ m = 2; // dimension of C
    ipc_ H_ne = 4; // number of elements of H
    ipc_ A_ne = 3; // number of elements of A
    ipc_ C_ne = 3; // number of elements of C
    ipc_ H_dense_ne = 6; // number of elements of H
    ipc_ A_dense_ne = 6; // number of elements of A
    ipc_ C_dense_ne = 3; // number of elements of C
    ipc_ H_row[] = {0, 1, 2, 2}; // row indices, NB lower triangle
    ipc_ H_col[] = {0, 1, 2, 0};
    ipc_ H_ptr[] = {0, 1, 2, 4};
    ipc_ A_row[] = {0, 0, 1};
    ipc_ A_col[] = {0, 1, 2};
    ipc_ A_ptr[] = {0, 2, 3};
    ipc_ C_row[] = {0, 1, 1}; // row indices, NB lower triangle
    ipc_ C_col[] = {0, 0, 1};
    ipc_ C_ptr[] = {0, 1, 3};
    rpc_ H_val[] = {1.0, 2.0, 3.0, 1.0};
    rpc_ A_val[] = {2.0, 1.0, 1.0};
    rpc_ C_val[] = {4.0, 1.0, 2.0};
    rpc_ H_dense[] = {1.0, 0.0, 2.0, 1.0, 0.0, 3.0};
    rpc_ A_dense[] = {2.0, 1.0, 0.0, 0.0, 0.0, 1.0};
    rpc_ C_dense[] = {4.0, 1.0, 2.0};
    rpc_ H_diag[] = {1.0, 1.0, 2.0};
    rpc_ C_diag[] = {4.0, 2.0};
    rpc_ H_scid[] = {2.0};
    rpc_ C_scid[] = {2.0};

    char st = ' ';
    ipc_ status;

    printf(" C sparse matrix indexing\n\n");

    printf(" basic tests of storage formats\n\n");

    for( ipc_ d=1; d <= 7; d++){

        // Initialize SBLS
        sbls_initialize( &data, &control, &status );
        //control.print_level = 1;
        //control.sls_control.print_level = 1;
        control.preconditioner = 2;
        control.factorization = 2;
        control.get_norm_residual = true;
        strcpy(control.symmetric_linear_solver, "sytr ") ;
        strcpy(control.definite_linear_solver, "sytr ") ;

        // Set user-defined control options
        control.f_indexing = false; // C sparse matrix indexing

```

--------------------------------

### DQP example calls: storage switch and shifted least-distance setup

Source: https://ralna.github.io/galahad_docs/html/C/dqp.html

The first block shows a switch over storage types (coordinate, sparse_by_rows, dense, diagonal, scaled_identity, identity) that imports and solves a quadratic program with dqp_import and dqp_solve_qp, then prints results and terminates. The second block initializes the shifted least-distance interface, sets Fortran indexing, and defines starting vectors and shifted data.

```c
switch(d){
            case 1: // sparse co-ordinate storage
                st = 'C';
                dqp_import( &control, &data, &status, n, m,
                           "coordinate", H_ne, H_row, H_col, NULL,
                           "coordinate", A_ne, A_row, A_col, NULL );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;
            printf(" case %1" d_ipc_ " break\n",d);
            case 2: // sparse by rows
                st = 'R';
                dqp_import( &control, &data, &status, n, m,
                            "sparse_by_rows", H_ne, NULL, H_col, H_ptr,
                            "sparse_by_rows", A_ne, NULL, A_col, A_ptr );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;
            case 3: // dense
                st = 'D';
                ipc_ H_dense_ne = 6; // number of elements of H
                ipc_ A_dense_ne = 6; // number of elements of A
                rpc_ H_dense[] = {1.0, 0.0, 1.0, 0.0, 0.0, 1.0};
                rpc_ A_dense[] = {2.0, 1.0, 0.0, 0.0, 1.0, 1.0};
                dqp_import( &control, &data, &status, n, m,
                            "dense", H_ne, NULL, NULL, NULL,
                            "dense", A_ne, NULL, NULL, NULL );
                dqp_solve_qp( &data, &status, n, m, H_dense_ne, H_dense, g, f,
                              A_dense_ne, A_dense, c_l, c_u, x_l, x_u,
                              x, c, y, z, x_stat, c_stat );
                break;
            case 4: // diagonal
                st = 'L';
                dqp_import( &control, &data, &status, n, m,
                            "diagonal", H_ne, NULL, NULL, NULL,
                            "sparse_by_rows", A_ne, NULL, A_col, A_ptr );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;

            case 5: // scaled identity
                st = 'S';
                dqp_import( &control, &data, &status, n, m,
                            "scaled_identity", H_ne, NULL, NULL, NULL,
                            "sparse_by_rows", A_ne, NULL, A_col, A_ptr );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;
            case 6: // identity
                st = 'I';
                dqp_import( &control, &data, &status, n, m,
                            "identity", H_ne, NULL, NULL, NULL,
                            "sparse_by_rows", A_ne, NULL, A_col, A_ptr );
                dqp_solve_qp( &data, &status, n, m, H_ne, H_val, g, f,
                              A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                              x_stat, c_stat );
                break;
            }
        dqp_information( &data, &inform, &status );

        if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_f.h
#include "galahad_pquad_f.h"
#else
            printf("%c:%6" d_ipc_ " iterations. Optimal objective "
                   "value = %.2f status = %1" d_ipc_ "\n",
                   st, inform.iter, inform.obj, inform.status);
#endif
        }else{
            printf("%c: DQP_solve exit status = %1" d_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
        dqp_terminate( &data, &control, &inform );
    }
```

```c
// test shifted least-distance interface
    for( ipc_ d=1; d <= 1; d++){

        // Initialize DQP
        dqp_initialize( &data, &control, &status );

        // Set user-defined control options
        control.f_indexing = true; // Fortran sparse matrix indexing

        // Start from 0
        rpc_ x[] = {0.0,0.0,0.0};
        rpc_ y[] = {0.0,0.0};
        rpc_ z[] = {0.0,0.0,0.0};

        // Set shifted least-distance data

        rpc_ w[] = {1.0,1.0,1.0};
        rpc_ x_0[] = {0.0,0.0,0.0};
```

--------------------------------

### Solving Quadratic Programs with QPA in Julia

Source: https://ralna.github.io/galahad_docs/html/Julia/qpa.html

This example demonstrates the setup and execution of the QPA solver using different matrix storage formats. It defines problem data, initializes the solver, and iterates through various storage configurations.

```julia
# test_qpa.jl
# Simple code to test the Julia interface to QPA

using GALAHAD
using Test
using Printf
using Accessors
using Quadmath

function test_qpa(::Type{T}, ::Type{INT}; sls::String="sytr") where {T,INT}
  # Derived types
  data = Ref{Ptr{Cvoid}}()
  control = Ref{qpa_control_type{T,INT}}()
  inform = Ref{qpa_inform_type{T,INT}}()

  # Set problem data
  n = INT(3)  # dimension
  m = INT(2)  # number of general constraints
  H_ne = INT(3)  # Hesssian elements
  H_row = INT[1, 2, 3]  # row indices, NB lower triangle
  H_col = INT[1, 2, 3]  # column indices, NB lower triangle
  H_ptr = INT[1, 2, 3, 4]  # row pointers
  H_val = T[1.0, 1.0, 1.0]  # values
  g = T[0.0, 2.0, 0.0]  # linear term in the objective
  f = one(T)  # constant term in the objective
  rho_g = T(0.1)  # penalty paramter for general constraints
  rho_b = T(0.1)  # penalty paramter for simple bound constraints
  A_ne = INT(4)  # Jacobian elements
  A_row = INT[1, 1, 2, 2]  # row indices
  A_col = INT[1, 2, 2, 3]  # row indices
  A_ptr = INT[1, 3, 5]  # row pointers
  A_val = T[2.0, 1.0, 1.0, 1.0]  # values
  c_l = T[1.0, 2.0]  # constraint lower bound
  c_u = T[2.0, 2.0]  # constraint upper bound
  x_l = T[-1.0, -Inf, -Inf]  # variable lower bound
  x_u = T[1.0, Inf, 2.0]  # variable upper bound

  # Set output storage
  c = zeros(T, m) # constraint values
  x_stat = zeros(INT, n) # variable status
  c_stat = zeros(INT, m) # constraint status
  st = ' '
  status = Ref{INT}()

  @printf(" Fortran sparse matrix indexing\n\n")
  @printf(" basic tests of qp storage formats\n\n")

  for d in 1:7
    # Initialize QPA
    qpa_initialize(T, INT, data, control, status)

    # Linear solvers
    @reset control[].symmetric_linear_solver = galahad_linear_solver(sls)

    # Start from 0
    x = T[0.0, 0.0, 0.0]
    y = T[0.0, 0.0]
    z = T[0.0, 0.0, 0.0]

    # sparse co-ordinate storage
    if d == 1
      st = 'C'
      qpa_import(T, INT, control, data, status, n, m,
                 "coordinate", H_ne, H_row, H_col, C_NULL,
                 "coordinate", A_ne, A_row, A_col, C_NULL)

      qpa_solve_qp(T, INT, data, status, n, m, H_ne, H_val, g, f,
                   A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                   x_stat, c_stat)
    end

    # sparse by rows
    if d == 2
      st = 'R'
      qpa_import(T, INT, control, data, status, n, m,
                 "sparse_by_rows", H_ne, C_NULL, H_col, H_ptr,
                 "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)

      qpa_solve_qp(T, INT, data, status, n, m, H_ne, H_val, g, f,
                   A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                   x_stat, c_stat)
    end

    # dense
    if d == 3
      st = 'D'
      H_dense_ne = 6 # number of elements of H
      A_dense_ne = 6 # number of elements of A
      H_dense = T[1.0, 0.0, 1.0, 0.0, 0.0, 1.0]
      A_dense = T[2.0, 1.0, 0.0, 0.0, 1.0, 1.0]

      qpa_import(T, INT, control, data, status, n, m,
                 "dense", H_ne, C_NULL, C_NULL, C_NULL,
                 "dense", A_ne, C_NULL, C_NULL, C_NULL)

      qpa_solve_qp(T, INT, data, status, n, m, H_dense_ne, H_dense, g, f,
                   A_dense_ne, A_dense, c_l, c_u, x_l, x_u,
                   x, c, y, z, x_stat, c_stat)
    end

    # diagonal
    if d == 4
      st = 'L'
      qpa_import(T, INT, control, data, status, n, m,
                 "diagonal", H_ne, C_NULL, C_NULL, C_NULL,
                 "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)

      qpa_solve_qp(T, INT, data, status, n, m, H_ne, H_val, g, f,
                   A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                   x_stat, c_stat)
    end

    # scaled identity
    if d == 5
      st = 'S'
      qpa_import(T, INT, control, data, status, n, m,
                 "scaled_identity", H_ne, C_NULL, C_NULL, C_NULL,
                 "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)

      qpa_solve_qp(T, INT, data, status, n, m, H_ne, H_val, g, f,
                   A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                   x_stat, c_stat)
    end

    # identity
    if d == 6
      st = 'I'
      qpa_import(T, INT, control, data, status, n, m,
                 "identity", H_ne, C_NULL, C_NULL, C_NULL,
                 "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)

      qpa_solve_qp(T, INT, data, status, n, m, H_ne, H_val, g, f,
                   A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
                   x_stat, c_stat)
    end

    # zero
    if d == 7
      st = 'Z'
      qpa_import(T, INT, control, data, status, n, m,
                 "zero", H_ne, C_NULL, C_NULL, C_NULL,
                 "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)
```

--------------------------------

### EXPO multi-precision C example: setup and callbacks

Source: https://ralna.github.io/galahad_docs/html/C/expo_multi.html

Shows the C setup for a multi-precision EXPO call using the _s_64 (single precision, 64-bit integer) interface. It sets problem dimensions, sparse Jacobian/Hessian indices, bounds, and control options, then imports, solves, prints status, and terminates. The second block defines the user-supplied fc, gj, and hl callbacks that compute the objective, constraints, gradient/Jacobian, and Hessian of the Lagrangian for the example problem.

```c
    // Set problem data
    int64_t n_64 = 2; // # variables
    int64_t m_64 = 5; // # constraints
    int64_t j_ne_64 = 10; // Jacobian elements
    int64_t h_ne_64 = 2; // Hesssian elements
    // 0-based indices
    int64_t J_row_64[] = {0, 0, 1, 1, 2, 2, 3, 3, 4, 4}; // Jacobian J
    int64_t J_col_64[] = {0, 1, 0, 1, 0, 1, 0, 1, 0, 1};
    int64_t H_row_64[] = {0, 1}; // Hessian H, NB lower triangle
    int64_t H_col_64[] = {0, 1};

    // Set storage
    float y_s[m]; // multipliers
    float z_s[n]; // dual variables
    float c_s[m]; // constraints
    float gl_s[n]; // gradients
    float xl_s[] = {-50.0, -50.0}; // lower variable bounds
    float xu_s[] = {50.0, 50.0}; // upper variable bounds
    float cl_s[] = {0.0, 0.0, 0.0, 0.0, 0.0}; // lower constraint bounds
    float cu_s[] = {INFINITY, INFINITY, INFINITY,
                   INFINITY, INFINITY}; // upper constraint bounds
    float x_s[] = {3.0,1.0}; // starting point
    int64_t status_64;

    // Initialize EXPO
    expo_initialize_s_64( &data_s_64, &control_s_64, &inform_s_64 );

    // Set user-defined control options
    control_s_64.f_indexing = false; // C sparse matrix indexing
    //control_s_64.print_level = 1;
    control_s_64.max_it = 20;
    control_s_64.max_eval = 100;
    control_s_64.stop_abs_p = 0.001;
    control_s_64.stop_abs_d = 0.001;
    control_s_64.stop_abs_c = 0.001;
    control_s_64.tru_control.error = 0;

    // Solve the problem
    expo_import_s_64( &control_s_64, &data_s_64, &status_64, n_64, m_64,
                 "coordinate", j_ne_64, J_row_64, J_col_64, NULL,
                 "coordinate", h_ne_64, H_row_64, H_col_64, NULL );
    expo_solve_hessian_direct_s_64( &data_s_64, &userdata_s_64, &status_64,
                                    n_64, m_64, j_ne_64, h_ne_64,
                                    cl_s, cu_s, xl_s, xu_s,
                                    x_s, y_s, z_s, c_s, gl_s,
                                    fc_s_64, gj_s_64, hl_s_64 );
    expo_information_s_64( &data_s_64, &inform_s_64, &status_64 );

    if(inform_s_64.status == 0){
      printf("%6li iterations. Optimal objective value = %.2f status = %1li\n",
          inform_s_64.iter, inform_s_64.obj, inform_s_64.status);
    }else{
      printf("EXPO_solve exit status = %1li\n", inform_s_64.status);
    }
    // Delete internal workspace
    expo_terminate_s_64( &data_s_64, &control_s_64, &inform_s_64 );

}
```

```c
// compute the function and constraint values
int fc( int32_t n, int32_t m, const double x[], double *f, double c[],
         const void *userdata ){
    struct userdata_type *myuserdata = (struct userdata_type *) userdata;
    double p = myuserdata->p;
    *f = pow(x[0],2.0) + pow(x[1],2.0);
    c[0] = x[0] + x[1] - 1.0;
    c[1] = pow(x[0],2.0) + pow(x[1],2.0) - 1.0;
    c[2] = p * pow(x[0],2.0) + pow(x[1],2.0) - p;
    c[3] = pow(x[0],2.0) - x[1];
    c[4] = pow(x[1],2.0) - x[0];
    return 0;
}

// compute the gradient and Jacobian
int gj( int32_t n, int32_t m, int32_t jne, const double x[], double g[],
        double jval[], const void *userdata ){
    struct userdata_type *myuserdata = (struct userdata_type *) userdata;
    double p = myuserdata->p;
    g[0] = 2.0 * x[0];
    g[1] = 2.0 * x[1];
    jval[0] = 1.0;
    jval[1] = 1.0;
    jval[2] = 2.0 * x[0];
    jval[3] = 2.0 * x[1];
    jval[4] = 2.0 * p * x[0];
    jval[5] = 2.0 * x[1];
    jval[6] = 2.0 * x[0];
    jval[7] = - 1.0;
    jval[8] = - 1.0;
    jval[9] = 2.0 * x[1];
    return 0;
}

// compute the Hessian of the Lagrangian
int hl( int32_t n, int32_t m, int32_t hne, const double x[], const double y[],
         double hval[], const void *userdata ){
    struct userdata_type *myuserdata = (struct userdata_type *) userdata;
    double p = myuserdata->p;
    hval[0] = 2.0 - 2.0 * (y[1] + p * y[2] + y[3]);
    hval[1] = 2.0 - 2.0 * (y[1] + y[2] + y[4]);
    return 0;
}
```

--------------------------------

### Solve a quadratic program with GALAHAD lsqp in Python

Source: https://ralna.github.io/galahad_docs/html/Python/lsqp.html

Complete example that initializes lsqp, loads problem data, solves a quadratic program, and retrieves solution information. Demonstrates setting options, providing starting values, and handling output status.

```python
from galahad import lsqp
import numpy as np
np.set_printoptions(precision=4,suppress=True,floatmode='fixed')
print("\n** python test: lsqp")

# set parameters
n = 3
m = 2
infinity = float("inf")

#  describe objective function

f = 1.0
g = np.array([0.0,2.0,0.0])
w = np.array([1.0,1.0,1.0])
x0 = np.array([1.0,1.0,1.0])

#  describe constraints

A_type = 'coordinate'
A_ne = 4
A_row = np.array([0,0,1,1])
A_col = np.array([0,1,1,2])
A_ptr = None
A_val = np.array([2.0,1.0,1.0,1.0])
c_l = np.array([1.0,2.0])
c_u = np.array([2.0,2.0])
x_l = np.array([-1.0,-infinity,-infinity])
x_u = np.array([1.0,infinity,2.0])

# allocate internal data and set default options
options = lsqp.initialize()

# set some non-default options
options['print_level'] = 0
#print("options:", options)

# load data (and optionally non-default options)
lsqp.load(n, m, A_type, A_ne, A_row, A_col, A_ptr, options)

#  provide starting values (not crucial)

x = np.array([0.0,0.0,0.0])
y = np.array([0.0,0.0])
z = np.array([0.0,0.0,0.0])

# find optimum of qp
#print("\nsolve lsqp")
x, c, y, z, x_stat, c_stat \
  = lsqp.solve_qp(n, m, f, g, w, x0, A_ne, A_val,
                 c_l, c_u, x_l, x_u, x, y, z)
print(" x:",x)
print(" c:",c)
print(" y:",y)
print(" z:",z)
print(" x_stat:",x_stat)
print(" c_stat:",c_stat)

# get information
inform = lsqp.information()
print(" f: %.4f" % inform['obj'])
print('** lsqp exit status:', inform['status'])

# deallocate internal data

lsqp.terminate()
```
Tokens:
Raw