Skip to content

Commit

Permalink
bug fix for ssids enquire
Browse files Browse the repository at this point in the history
  • Loading branch information
@dalekopera
dalekopera committed Jan 24, 2025
1 parent 5230165 commit 013c779
Show file tree
Hide file tree
Showing 9 changed files with 74 additions and 69 deletions.
23 changes: 16 additions & 7 deletions include/ssids_cpu_NumericSubtree.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
* \copyright 2016 The Science and Technology Facilities Council (STFC)
* \licence BSD licence, see LICENCE file for details
* \author Jonathan Hogg
* \version GALAHAD 5.2 - 2025-01-22 AT 12:00 GMT
* \version GALAHAD 5.2 - 2025-01-24 AT 14:00 GMT
*/

#pragma once
Expand Down Expand Up @@ -485,7 +485,8 @@ public:
/* 1x1 pivot */
if(piv_order) {
piv_order[nodes_[ni].perm[i]-1] = (piv++);
//printf(" 1x1 pivot order %d = %d\n", nodes_[ni].perm[i]-1, piv_order[nodes_[ni].perm[i]-1]);
// printf(" 1x1 pivot order %d = %d\n", nodes_[ni].perm[i]-1,
// piv_order[nodes_[ni].perm[i]-1]);
}
if(d) {
// printf("in = %i d(1,1) = %.1f\n", 2*i+0, dptr[2*i+0]);
Expand All @@ -498,15 +499,17 @@ public:
/* 2x2 pivot */
if(piv_order) {
piv_order[nodes_[ni].perm[i]-1] = -(piv++);
//printf(" 2x2 pivot order %d = %d\n", nodes_[ni].perm[i]-1, piv_order[nodes_[ni].perm[i]-1]);
// printf(" 2x2 pivot order %d = %d\n", nodes_[ni].perm[i]-1,
// piv_order[nodes_[ni].perm[i]-1]);
piv_order[nodes_[ni].perm[i+1]-1] = -(piv++);
//printf(" 2x2 pivot order %d = %d\n", nodes_[ni].perm[i+1]-1, piv_order[nodes_[ni].perm[i+1]-1]);
// printf(" 2x2 pivot order %d = %d\n", nodes_[ni].perm[i+1]-1,
// piv_order[nodes_[ni].perm[i+1]-1]);
}
if(d) {
// printf("in = %i d(1,1) = %.1f\n", 2*i+0, dptr[2*i+0]);
// printf("in = %i d(2,1) = %.1f\n", 2*i+1, dptr[2*i+1]);
// printf("in = %i d(1,2) = %.1f\n", 2*i+2, dptr[2*i+2]);
/// printf("in = %i d(2,2) = %.1f\n", 2*i+3, dptr[2*i+3]);
// printf("in = %i d(2,2) = %.1f\n", 2*i+3, dptr[2*i+3]);
*(d++) = dptr[2*i+0];
*(d++) = dptr[2*i+1];
*(d++) = dptr[2*i+3]; /* not 2*i+2 as stated ?? */
Expand All @@ -515,7 +518,13 @@ public:
i+=2;
}
}

}
// printf("piv_order: ");
// for(ipc_ i=0; i<4; i++) {
// printf(" %d", piv_order[i]);
// }
// printf("\n");
}
}

Expand All @@ -536,15 +545,15 @@ public:
#endif
/* 1x1 pivot */
dptr[2*i+0] = *(d++);
d++; /* correct incremenet */
d++; /* bug fix - correct increment */
// printf("in = %i d(1,1) = %.1f\n", 2*i+0, dptr[2*i+0]);
i+=1;
} else {
/* 2x2 pivot */
dptr[2*i+0] = *(d++);
dptr[2*i+1] = *(d++);
dptr[2*i+3] = *(d++);
d++; /* correct incremenet */
d++; /* bug fix - correct increment */
// printf("in = %i d(1,1) = %.1f\n", 2*i+0, dptr[2*i+0]);
// printf("in = %i d(2,1) = %.1f\n", 2*i+1, dptr[2*i+1]);
// printf("in = %i d(1,2) = %.1f\n", 2*i+2, dptr[2*i+2]);
Expand Down
9 changes: 8 additions & 1 deletion src/cqp/C/cqpt.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_cqp.h"
Expand Down Expand Up @@ -51,11 +52,14 @@ int main(void) {
for( ipc_ d=1; d <= 7; d++){

// Initialize CQP
printf("initialize\n");
cqp_initialize( &data, &control, &status );

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

control.print_level = 3;
// strcpy(control.fdc_control.symmetric_linear_solver, "sytr ");
control.fdc_control.print_level = 3;
// Start from 0
rpc_ x[] = {0.0,0.0,0.0};
rpc_ y[] = {0.0,0.0};
Expand All @@ -64,9 +68,11 @@ int main(void) {
switch(d){
case 1: // sparse co-ordinate storage
st = 'C';
printf("import\n");
cqp_import( &control, &data, &status, n, m,
"coordinate", H_ne, H_row, H_col, NULL,
"coordinate", A_ne, A_row, A_col, NULL );
printf("solve\n");
cqp_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 );
Expand Down Expand Up @@ -135,6 +141,7 @@ int main(void) {


}
printf("info\n");
cqp_information( &data, &inform, &status );

if(inform.status == 0){
Expand Down
4 changes: 3 additions & 1 deletion src/cqp/cqpti.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,9 @@ PROGRAM GALAHAD_CQP_interface_test

DO data_storage_type = 1, 7
CALL CQP_initialize( data, control, inform )
CALL WHICH_sls( control )
control%print_level = 3
control%fdc_control%print_level = 4
! CALL WHICH_sls( control )
X = 0.0_rp_ ; Y = 0.0_rp_ ; Z = 0.0_rp_ ! start from zero
SELECT CASE ( data_storage_type )
CASE ( 1 ) ! sparse co-ordinate storage
Expand Down
10 changes: 7 additions & 3 deletions src/eqp/eqpdt.output
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,10 @@

basic tests of qp storage formats

At line 1306 of file eqp.F90
Fortran runtime error: Index '0' of dimension 1 of array 'data' outside of expected range (1:2)
2 0
C: 1 cg iterations. Optimal objective value = 0.44 status = 0
R: 1 cg iterations. Optimal objective value = 0.44 status = 0
D: 1 cg iterations. Optimal objective value = 0.44 status = 0
L: 1 cg iterations. Optimal objective value = 0.44 status = 0
S: 1 cg iterations. Optimal objective value = 0.44 status = 0
I: 1 cg iterations. Optimal objective value = 0.44 status = 0
W: 1 cg iterations. Optimal objective value = 0.44 status = 0
17 changes: 10 additions & 7 deletions src/fdc/fdc.F90
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
! THIS VERSION: GALAHAD 5.2 - 2025-01-22 AT 14:10 GMT.
! THIS VERSION: GALAHAD 5.2 - 2025-01-24 AT 10:00 GMT.

#include "galahad_modules.h"

Expand Down Expand Up @@ -967,12 +967,15 @@ SUBROUTINE FDC_find_dependent_sls( n, m, A_val, A_col, A_ptr, C, &
RETURN
END IF

! write(6,"( ' solver is is ', A )" ) TRIM( inform%SLS_inform%solver )
! write(6,"( ' fdc pivots ', 10I6, /, ( 10I6 ) )" ) data%P( : data%K%n )
! write(6,"( ' fdc D_diag ', 5ES12.4, /, ( 12X, 5ES12.4 ) )" ) &
! data%D( 1, : data%K%n )
! write(6,"( ' fdc D_off ', 5ES12.4, /, ( 12X, 5ES12.4 ) )" ) &
! data%D( 2, : data%K%n )
IF ( out > 0 .AND. data%control%print_level >= 4 ) THEN
WRITE( out, "( ' solver is is ', A )" ) TRIM( inform%SLS_inform%solver )
WRITE( out, "( ' fdc pivots ', 10I6, /, ( 10I6 ) )" ) &
data%P( : data%K%n )
WRITE( out, "( ' fdc D_diag ', 5ES12.4, /, ( 12X, 5ES12.4 ) )" ) &
data%D( 1, : data%K%n )
WRITE( out, "( ' fdc D_off ', 5ES12.4, /, ( 12X, 5ES12.4 ) )" ) &
data%D( 2, : data%K%n )
END IF

! Compute the smallest and largest eigenvalues of the block diagonal factor

Expand Down
6 changes: 3 additions & 3 deletions src/fdc/fdct.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,16 +22,16 @@ PROGRAM GALAHAD_FDC_test !! to be expanded
control%use_sls = .TRUE.
control%symmetric_linear_solver = 'sytr'
! control%symmetric_linear_solver = 'ma57'
! control%symmetric_linear_solver = 'ssids'
control%symmetric_linear_solver = 'ssids'
CALL FDC_find_dependent( n, m, A_val, A_col, A_ptr, B, n_depen, DEPEN, &
data, control, inform ) ! Check for dependencies
WRITE( 6, "( ' linear solver used: ', A )" ) inform%SLS_inform%solver
IF ( inform%status == 0 ) THEN ! Successful return
IF ( n_depen == 0 ) THEN
WRITE( 6, "( ' FDC_find_dependent - no dependencies ' )" )
ELSE
WRITE( 6, "( ' FDC_find_dependent - dependent constraint(s):', 3I3 )") &
DEPEN( : n_depen )
WRITE( 6, "( ' FDC_find_dependent -', &
& ' indices of dependent constraint(s):', 3I3 )") DEPEN( : n_depen )
END IF
ELSE ! Error returns
WRITE( 6, "( ' FDC_find_dependent exit status = ', I6 ) " ) inform%status
Expand Down
37 changes: 0 additions & 37 deletions src/glrt/glrtdt.output
View file
Original file line number Diff line number Diff line change
@@ -1,41 +1,4 @@
MR = 00
MR = 00 glrt_solve_problem exit status = 0, f = -7.57
MR = 01
MR = 01 glrt_solve_problem exit status = 0, f = -3.34
MR = 10
MR = 10 glrt_solve_problem exit status = 0, f = -9.63
MR = 11

stopping tolerance = 1.0537E-07, sigma = 1.0000E+00

Iter objective pgnorm lambda gamma it info
0 0.00000000E+00 7.07E+00 - - - -
1 -7.56564165E+00 1.29E-01 1.84420637007342E+00 7.00E-02 2 0d
2 -7.56857169E+00 2.27E-02 1.84514276637603E+00 5.00E-01 3 0d
3 -7.56866504E+00 4.11E-03 1.84519734730030E+00 5.00E-01 2 0d
4 -7.56866811E+00 7.47E-04 1.84519998470893E+00 5.00E-01 2 0d
5 -7.56866821E+00 1.35E-04 1.84520009936154E+00 5.00E-01 2 0d
6 -7.56866821E+00 2.46E-05 1.84520010405336E+00 5.00E-01 2 0d
7 -7.56866821E+00 4.46E-06 1.84520010423805E+00 5.00E-01 2 0d
8 -7.56866821E+00 8.10E-07 1.84520010424513E+00 5.00E-01 2 0d
9 -7.56866821E+00 1.47E-07 1.84520010424513E+00 5.00E-01 1 0d
10 -7.56866821E+00 2.67E-08 1.84520010424513E+00 5.00E-01 1 0d

stopping tolerance = 1.4901E-07, sigma = 1.0000E+00

Iter objective pgnorm lambda gamma it info
0 0.00000000E+00 1.00E+01 - - - -
1 -9.61921723E+00 2.44E-01 1.74632172062310E+00 1.40E-01 2 0d
2 -9.62719764E+00 6.53E-02 1.74847790614645E+00 1.00E+00 3 0d
3 -9.62780907E+00 1.87E-02 1.74879664044715E+00 1.00E+00 3 0d
4 -9.62785977E+00 5.41E-03 1.74883686433044E+00 1.00E+00 2 0d
5 -9.62786401E+00 1.56E-03 1.74884138743299E+00 1.00E+00 2 0d
6 -9.62786436E+00 4.52E-04 1.74884186268817E+00 1.00E+00 2 0d
7 -9.62786439E+00 1.31E-04 1.74884191052577E+00 1.00E+00 2 0d
8 -9.62786439E+00 3.77E-05 1.74884191520131E+00 1.00E+00 2 0d
9 -9.62786439E+00 1.09E-05 1.74884191564860E+00 1.00E+00 2 0d
10 -9.62786439E+00 3.15E-06 1.74884191569070E+00 1.00E+00 2 0d
11 -9.62786439E+00 9.11E-07 1.74884191569462E+00 1.00E+00 2 0d
12 -9.62786439E+00 2.63E-07 1.74884191569462E+00 1.00E+00 1 0d
13 -9.62786439E+00 7.61E-08 1.74884191569462E+00 1.00E+00 1 0d
MR = 11 glrt_solve_problem exit status = 0, f = -5.02
1 change: 1 addition & 0 deletions src/ssids/fkeep.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -413,6 +413,7 @@ subroutine enquire_indef_cpu(akeep, fkeep, inform, piv_order, d)
if(present(piv_order)) then
do i = 1, akeep%n
piv_order( akeep%invp(i) ) = po(i)
! piv_order( i ) = po(i)
end do
endif

Expand Down
36 changes: 26 additions & 10 deletions src/ssids/ssids.F90
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
! THIS VERSION: GALAHAD 5.2 - 2025-01-22 AT 14:00 GMT.
! THIS VERSION: GALAHAD 5.2 - 2025-01-24 AT 14:00 GMT.
! (consistent with SPRAL up to issue #150)

#include "spral_procedures.h"
Expand Down Expand Up @@ -1361,17 +1361,33 @@ subroutine ssids_enquire_indef_precision(akeep, fkeep, options, inform, &
call fkeep%enquire_indef(akeep, inform, piv_order, d)
! bug fix to give 1-based indices
IF ( PRESENT( piv_order ) ) THEN
!write(6,"(' ssids: piv_order ', 7I6 )" ) piv_order( : akeep%n )
DO i = 1, akeep%n
po = piv_order( i )
! write(6,"(' ssids: piv_order ', 7I6 )" ) piv_order( : akeep%n )

! bug fix to determine what a C 0 index means
IF ( MOD( COUNT( piv_order( : akeep%n ) < 0 ), 2 ) == 0 ) THEN
DO i = 1, akeep%n
po = piv_order( i )
IF ( po >= 0 ) THEN
! bug fix as 0-based 0 is positive!
! IF ( po > 0 ) THEN
piv_order( i ) = po + 1
ELSE
piv_order( i ) = po - 1
END IF
END DO
ELSE
DO i = 1, akeep%n
po = piv_order( i )
! IF ( po >= 0 ) THEN
! bug fix as 0-based 0 is positive!
IF ( po > 0 ) THEN
piv_order( i ) = po + 1
ELSE
piv_order( i ) = po - 1
END IF
END DO
IF ( po > 0 ) THEN
piv_order( i ) = po + 1
ELSE
piv_order( i ) = po - 1
END IF
END DO
END IF
! write(6,"(' ssids: revised piv_order ', 7I6 )" ) piv_order( : akeep%n )
END IF
call inform%print_flag(options, context)
end subroutine ssids_enquire_indef_precision
Expand Down

0 comments on commit 013c779

Please sign in to comment.