--- Begin Message ---
Package: release.debian.org
Severity: normal
User: release.debian.org@packages.debian.org
Usertags: unblock
Please unblock package graywolf
Needed to fix release critical license violation. (#856705)
Basically two files (with license violation) are removed from upstream source.
Then the code is adapted to use libgsl instead.
Complete debdiff attached.
unblock graywolf/0.1.4+20170306gitecee764-1
-- System Information:
Debian Release: 8.7
APT prefers stable-updates
APT policy: (500, 'stable-updates'), (500, 'stable')
Architecture: amd64 (x86_64)
Kernel: Linux 3.16.0-4-amd64 (SMP w/1 CPU core)
Locale: LANG=en_US.UTF-8, LC_CTYPE=en_US.UTF-8 (charmap=UTF-8)
Shell: /bin/sh linked to /bin/dash
Init: systemd (via /run/systemd/system)
diff -Nru graywolf-0.1.4/CMakeLists.txt graywolf-0.1.4+20170306gitecee764/CMakeLists.txt
--- graywolf-0.1.4/CMakeLists.txt 2016-10-06 22:10:19.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/CMakeLists.txt 2017-03-06 20:47:22.000000000 +0100
@@ -1,8 +1,16 @@
cmake_minimum_required (VERSION 2.6)
project (graywolf)
+find_package(PkgConfig)
INCLUDE(CheckIncludeFiles)
+
+pkg_check_modules(GSL gsl)
+
+if (NOT GSL_FOUND)
+ MESSAGE(FATAL_ERROR "The development files for the GNU Scientific Library (libgsl) are required to build graywolf.")
+endif()
+
# Include RPATH in build so that ldconfig is not necessary after install
SET(CMAKE_SKIP_BUILD_RPATH FALSE)
SET(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
diff -Nru graywolf-0.1.4/debian/changelog graywolf-0.1.4+20170306gitecee764/debian/changelog
--- graywolf-0.1.4/debian/changelog 2016-10-06 22:20:06.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/debian/changelog 2017-03-06 20:48:51.000000000 +0100
@@ -1,3 +1,13 @@
+graywolf (0.1.4+20170306gitecee764-1) unstable; urgency=medium
+
+ * New upstream GIT HEAD with fixed license violation issue (Closes: #856705)
+ * debian/control:
+ - New build-dependencies: libgsl-dev, pkg-config
+ * Updated debian/libycadgraywolf0.symbols - removed the symbols in license
+ violating code.
+
+ -- Ruben Undheim <ruben.undheim@gmail.com> Mon, 06 Mar 2017 19:48:51 +0000
+
graywolf (0.1.4-1) unstable; urgency=low
* New upstream release
diff -Nru graywolf-0.1.4/debian/control graywolf-0.1.4+20170306gitecee764/debian/control
--- graywolf-0.1.4/debian/control 2016-10-06 22:20:06.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/debian/control 2017-03-06 20:48:51.000000000 +0100
@@ -5,7 +5,9 @@
Priority: optional
Build-Depends: debhelper (>= 9),
cmake,
- libx11-dev
+ pkg-config,
+ libx11-dev,
+ libgsl-dev
Standards-Version: 3.9.8
Vcs-Browser: https://anonscm.debian.org/cgit/debian-science/packages/graywolf.git
Vcs-Git: https://anonscm.debian.org/git/debian-science/packages/graywolf.git
diff -Nru graywolf-0.1.4/debian/libycadgraywolf0.symbols graywolf-0.1.4+20170306gitecee764/debian/libycadgraywolf0.symbols
--- graywolf-0.1.4/debian/libycadgraywolf0.symbols 2016-10-06 22:20:06.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/debian/libycadgraywolf0.symbols 2017-03-06 20:48:51.000000000 +0100
@@ -267,17 +267,6 @@
Ylist_sort@Base 0.1.1
Ylog_msg@Base 0.1.1
Ylog_start@Base 0.1.1
- Ymatrix_cofactors@Base 0.1.1
- Ymatrix_copy@Base 0.1.1
- Ymatrix_create@Base 0.1.1
- Ymatrix_disp@Base 0.1.1
- Ymatrix_eye@Base 0.1.1
- Ymatrix_free@Base 0.1.1
- Ymatrix_linv@Base 0.1.1
- Ymatrix_mult@Base 0.1.1
- Ymatrix_sub@Base 0.1.1
- Ymatrix_transpose@Base 0.1.1
- Ymatrix_zero@Base 0.1.1
Ymessage_close@Base 0.1.1
Ymessage_error_count@Base 0.1.1
Ymessage_flush@Base 0.1.1
@@ -389,8 +378,6 @@
Ystat_var@Base 0.1.1
Ystrclone@Base 0.1.1
Ystrparser@Base 0.1.1
- Ysvd_decompose@Base 0.1.1
- Ysvd_solve@Base 0.1.1
Ysystem@Base 0.1.1
Ytimer_elapsed@Base 0.1.1
Ytimer_start@Base 0.1.1
diff -Nru graywolf-0.1.4/debian/patches/01_no_rpath_in_debian_package.patch graywolf-0.1.4+20170306gitecee764/debian/patches/01_no_rpath_in_debian_package.patch
--- graywolf-0.1.4/debian/patches/01_no_rpath_in_debian_package.patch 2016-10-06 22:20:06.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/debian/patches/01_no_rpath_in_debian_package.patch 2017-03-06 20:48:51.000000000 +0100
@@ -8,13 +8,13 @@
CMakeLists.txt | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
-diff --git a/CMakeLists.txt b/CMakeLists.txt
-index 10e2b22..b8b6444 100644
---- a/CMakeLists.txt
-+++ b/CMakeLists.txt
-@@ -3,11 +3,11 @@ project (graywolf)
-
- INCLUDE(CheckIncludeFiles)
+Index: graywolf/CMakeLists.txt
+===================================================================
+--- graywolf.orig/CMakeLists.txt
++++ graywolf/CMakeLists.txt
+@@ -11,11 +11,11 @@ if (NOT GSL_FOUND)
+ MESSAGE(FATAL_ERROR "The development files for the GNU Scientific Library (libgsl) are required to build graywolf.")
+ endif()
-# Include RPATH in build so that ldconfig is not necessary after install
-SET(CMAKE_SKIP_BUILD_RPATH FALSE)
diff -Nru graywolf-0.1.4/debian/patches/02_dont_create_link.patch graywolf-0.1.4+20170306gitecee764/debian/patches/02_dont_create_link.patch
--- graywolf-0.1.4/debian/patches/02_dont_create_link.patch 2016-10-06 22:20:06.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/debian/patches/02_dont_create_link.patch 2017-03-06 20:48:51.000000000 +0100
@@ -8,11 +8,11 @@
CMakeLists.txt | 1 -
1 file changed, 1 deletion(-)
-diff --git a/CMakeLists.txt b/CMakeLists.txt
-index b8b6444..2411c46 100644
---- a/CMakeLists.txt
-+++ b/CMakeLists.txt
-@@ -14,7 +14,6 @@ CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DIR}/config.h.in ${CMAKE_CURRENT_BINARY_DI
+Index: graywolf/CMakeLists.txt
+===================================================================
+--- graywolf.orig/CMakeLists.txt
++++ graywolf/CMakeLists.txt
+@@ -22,7 +22,6 @@ CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DI
install(DIRECTORY flow DESTINATION lib/graywolf/bin)
diff -Nru graywolf-0.1.4/src/twmc/CMakeLists.txt graywolf-0.1.4+20170306gitecee764/src/twmc/CMakeLists.txt
--- graywolf-0.1.4/src/twmc/CMakeLists.txt 2016-10-06 22:10:19.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/src/twmc/CMakeLists.txt 2017-03-06 20:47:22.000000000 +0100
@@ -5,6 +5,7 @@
target_link_libraries(TimberWolfMC ${CMAKE_BINARY_DIR}/src/Ylib/libycadgraywolf.so)
target_link_libraries(TimberWolfMC X11)
target_link_libraries(TimberWolfMC m)
+target_link_libraries(TimberWolfMC ${GSL_LIBRARIES})
INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/include ${CMAKE_BINARY_DIR}/include .)
diff -Nru graywolf-0.1.4/src/twmc/wireratio.c graywolf-0.1.4+20170306gitecee764/src/twmc/wireratio.c
--- graywolf-0.1.4/src/twmc/wireratio.c 2016-10-06 22:10:19.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/src/twmc/wireratio.c 2017-03-06 20:47:22.000000000 +0100
@@ -60,10 +60,39 @@
#include <dens.h>
#include <yalecad/debug.h>
#include <yalecad/file.h>
-#include <yalecad/linalg.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_linalg.h>
+
+gsl_matrix_disp( mptr, rows, cols )
+gsl_matrix *mptr ;
+int rows, cols;
+{
+ INT i, j ;
+
+ for( i=0; i < rows; i++ ){
+ for( j=0; j < cols; j++ ){
+ fprintf( stderr, "% 4.4le ", gsl_matrix_get(mptr, i, j)) ;
+ }
+ fprintf( stderr, "\n" ) ;
+ }
+ fprintf( stderr, "\n" ) ;
+} /* end gsl_matrix_disp */
+
+gsl_vector_disp( vptr, rows )
+gsl_vector *vptr ;
+int rows;
+{
+ INT i;
+
+ for( i=0; i < rows; i++ ){
+ fprintf( stderr, "% 4.4le ", gsl_vector_get(vptr, i)) ;
+ }
+ fprintf( stderr, "\n" ) ;
+} /* end gsl_vector_disp */
static set_pins( A, center, loc, tile_side, sidepins, count )
-YMPTR A ; /* the matrix holding x y positions */
+gsl_matrix *A ; /* the matrix holding x y positions */
INT center ;
INT loc ;
INT tile_side ;
@@ -75,15 +104,15 @@
if( sidepins ){
side = find_tile_side( center, loc, tile_side ) ;
if( side ){
- YMATRIX( A, count, 6 ) = (DOUBLE) sidepins[side] ;
+ gsl_matrix_set(A, count, 5, (DOUBLE) sidepins[side]);
} else {
M( ERRMSG, "adapt_wire_estimator",
"Trouble finding pinside - defaulting to 0.\n" ) ;
- YMATRIX( A, count, 6 ) = 0.0 ;
+ gsl_matrix_set(A, count, 5, 0.0);
}
} else {
/* no pins for cell */
- YMATRIX( A, count, 6 ) = 0.0 ;
+ gsl_matrix_set(A, count, 5, 0.0);
}
} /* end set_pins */
@@ -95,12 +124,13 @@
INT xc, yc ; /* cell center */
INT *sidepins ; /* array holding #pins for side */
INT l, r, b, t ; /* the global position of the rtiles */
+ INT solved ; /* status of gsl_linalg_SV_solve */
INT *find_pin_sides() ; /* find number of pins on all sides */
char filename[LRECL] ; /* output the results of the SVD fit */
FILE *fp ; /* write out the results */
- YMPTR A ; /* the matrix holding x y positions */
- YMPTR B ; /* the resulting global routing space*/
- YMPTR Xret ; /* the solution to At*A*x = At*B */
+ gsl_matrix *A ; /* the matrix holding x y positions */
+ gsl_vector *B ; /* the resulting global routing space*/
+ gsl_vector *Xret ; /* the solution to At*A*x = At*B */
DOUBLE x, y ; /* cell placement */
DOUBLE lf, rf, bf, tf ; /* the global position of the rtiles */
DOUBLE xlength, ylength;/* length of side */
@@ -108,6 +138,11 @@
CELLBOXPTR cptr ; /* current pointer to cell */
RTILEBOXPTR rptr ; /* traverse tiles */
+ gsl_matrix *U ;
+ gsl_matrix *V ;
+ gsl_vector *S ;
+ gsl_vector *work ;
+
if(!(routingTilesG)){
return ;
}
@@ -124,12 +159,12 @@
}
if( count < 6 ){
/* now we can size the matrices */
- A = Ymatrix_create( 6, 6 ) ;
- B = Ymatrix_create( 6, 1 ) ;
+ A = gsl_matrix_alloc(6, 6);
+ B = gsl_vector_alloc(6);
} else {
/* now we can size the matrices */
- A = Ymatrix_create( count, 6 ) ;
- B = Ymatrix_create( count, 1 ) ;
+ A = gsl_matrix_alloc(count, 6);
+ B = gsl_vector_alloc(count);
}
/* initialize the pin counting routines */
@@ -151,8 +186,7 @@
/* switchbox densities are not accurate discard */
continue ;
}
- count++ ;
- YMATRIX( A, count, 1 ) = 1.0 ; /* for finding const */
+ gsl_matrix_set(A, count, 0, 1.0); /* for finding const */
l = xc + rptr->x1 ;
r = xc + rptr->x2 ;
b = yc + rptr->y1 ;
@@ -171,46 +205,47 @@
x = rf ;
y = (bf + tf) / 2.0 ;
set_pins( A, (b+t)/2, r, TILEL, sidepins, count ) ;
- YMATRIX( A, count, 2 ) = x ;
- YMATRIX( A, count, 3 ) = x * x ;
- YMATRIX( A, count, 4 ) = y ;
- YMATRIX( A, count, 5 ) = y * y ;
- YMATRIX( B, count, 1 ) = xrouting ;
+ gsl_matrix_set(A, count, 1, x);
+ gsl_matrix_set(A, count, 2, x * x);
+ gsl_matrix_set(A, count, 3, y);
+ gsl_matrix_set(A, count, 4, y * y);
+ gsl_vector_set(B, count, xrouting);
break ;
case TILET:
/* calculate x and y */
x = (lf + rf) / 2.0 ;
y = bf ;
set_pins( A, (l+r)/2, b, TILET, sidepins, count ) ;
- YMATRIX( A, count, 2 ) = x ;
- YMATRIX( A, count, 3 ) = x * x ;
- YMATRIX( A, count, 4 ) = y ;
- YMATRIX( A, count, 5 ) = y * y ;
- YMATRIX( B, count, 1 ) = yrouting ;
+ gsl_matrix_set(A, count, 1, x);
+ gsl_matrix_set(A, count, 2, x * x);
+ gsl_matrix_set(A, count, 3, y);
+ gsl_matrix_set(A, count, 4, y * y);
+ gsl_vector_set(B, count, yrouting);
break ;
case TILER:
/* calculate x and y */
x = lf ;
y = (bf + tf) / 2.0 ;
set_pins( A, (b+t)/2, l, TILER, sidepins, count ) ;
- YMATRIX( A, count, 2 ) = x ;
- YMATRIX( A, count, 3 ) = x * x ;
- YMATRIX( A, count, 4 ) = y ;
- YMATRIX( A, count, 5 ) = y * y ;
- YMATRIX( B, count, 1 ) = xrouting ;
+ gsl_matrix_set(A, count, 1, x);
+ gsl_matrix_set(A, count, 2, x * x);
+ gsl_matrix_set(A, count, 3, y);
+ gsl_matrix_set(A, count, 4, y * y);
+ gsl_vector_set(B, count, xrouting);
break ;
case TILEB:
/* calculate x and y */
x = (lf + rf) / 2.0 ;
y = tf ;
set_pins( A, (l+r)/2, t, TILEB, sidepins, count ) ;
- YMATRIX( A, count, 2 ) = x ;
- YMATRIX( A, count, 3 ) = x * x ;
- YMATRIX( A, count, 4 ) = y ;
- YMATRIX( A, count, 5 ) = y * y ;
- YMATRIX( B, count, 1 ) = yrouting ;
+ gsl_matrix_set(A, count, 1, x);
+ gsl_matrix_set(A, count, 2, x * x);
+ gsl_matrix_set(A, count, 3, y);
+ gsl_matrix_set(A, count, 4, y * y);
+ gsl_vector_set(B, count, yrouting);
break ;
} /* end switch */
+ count++ ;
}
if( sidepins ){
Yvector_free( sidepins, 1, sizeof(INT) ) ;
@@ -219,68 +254,98 @@
} /* end loop on cells */
/* now make sure we have at least 6 rows */
- if( count <= 1 ){
+ if( count < 1 ){
M( ERRMSG, "adapt_wire_estimator",
"Too few cells for TimberWolfMC. Must abort\n" ) ;
YexitPgm(PGMFAIL) ;
} else if( count < 6 ){
/* pad with zeros */
- for( ++count; count <= 6; count++ ){
- for( i = 1; i <= 6; i++ ){
- YMATRIX( A, count, i ) = YMATRIX( A, 1, i ) ;
+ for( ; count < 6; count++ ){
+ for( i = 0; i < 6; i++ ){
+ gsl_matrix_set(A, count, i, gsl_matrix_get(A, 0, i));
}
- YMATRIX( B, count, 1 ) = YMATRIX( B, 1, 1 ) ;
+ gsl_vector_set(B, count, gsl_vector_get(B, 0));
}
}
D( "TWMC/awe/init",
fprintf( stderr,"\nA:\n" ) ;
- Ymatrix_disp( A ) ;
+ gsl_matrix_disp( A, count < 6 ? 6 : count, 6 ) ;
fprintf( stderr,"\nB:\n" ) ;
- Ymatrix_disp( B ) ;
+ gsl_vector_disp( B, count < 6 ? 6 : count ) ;
) ;
/* now solve using SVD */
- if( Ysvd_solve( A, B, &Xret )){
+ U = gsl_matrix_alloc((count < 6) ? 6 : count, 6);
+ V = gsl_matrix_alloc(6, 6);
+ S = gsl_vector_alloc(6);
+ work = gsl_vector_alloc(6);
+ gsl_matrix_memcpy(U, A);
+ gsl_linalg_SV_decomp(U, V, S, work);
+ solved = gsl_linalg_SV_solve(U, V, S, B, Xret);
+
+ if( solved ){
sprintf( filename, "%s.mest", cktNameG ) ;
fp = TWOPEN( filename, "w", ABORT ) ;
fprintf( fpoG, "\nThe results of the SVD fit are:\n" ) ;
- for( i = 1; i <= 6; i++ ){
- HPO( fp, YMATRIX(Xret,i,1) ) ;
- HPO( fpoG, YMATRIX(Xret,i,1) ) ;
+ for( i = 0; i < 6; i++ ){
+ HPO( fp, gsl_vector_get(Xret, i));
+ HPO( fpoG, gsl_vector_get(Xret, i));
}
TWCLOSE( fp ) ;
}
+ gsl_matrix_free(U);
+ gsl_matrix_free(V);
+ gsl_vector_free(S);
+ gsl_vector_free(work);
D( "TWMC/awe/answer",
- YMPTR AX ; /* multiply the answer */
- YMPTR R ;
- YMPTR Q ;
+ gsl_vector *AX ; /* multiply the answer */
+ gsl_vector *R ;
+ gsl_matrix *Q ;
INT c ;
+ DOUBLE d ;
+ DOUBLE dx ;
+
+ AX = gsl_vector_alloc(count < 6 ? count : 6);
+ R = gsl_vector_alloc(count < 6 ? count : 6);
+
+ /* Compute AX = A * Xret */
+ for ( c = 0; c < (count < 6) ? 6 : count; c++ ) {
+ d = 0.0;
+ dx = gsl_vector_get(Xret, c);
+ for ( i = 0; i < 6; i++ ) {
+ d += gsl_matrix_get(A, c, i) * dx;
+ gsl_set_vector(AX, c, d);
+ }
- /* now output the differences */
- AX = Ymatrix_mult( A, Xret ) ;
- R = Ymatrix_sub( AX, B ) ;
- Q = Ymatrix_create( count, 7 ) ;
- for( c = 1; c <= count; c++ ){
- YMATRIX( Q, c, 1 ) = YMATRIX( A, c, 2 ) ; /* x */
- YMATRIX( Q, c, 2 ) = YMATRIX( A, c, 4 ) ; /* y */
- YMATRIX( Q, c, 3 ) = YMATRIX( A, c, 6 ) ; /* pins */
- YMATRIX( Q, c, 4 ) = YMATRIX( B, c, 1 ) ; /* B */
- YMATRIX( Q, c, 5 ) = YMATRIX( AX, c, 1 ) ; /* AX */
- YMATRIX( Q, c, 6 ) = YMATRIX( R, c, 1 ) ; /* B - AX */
- YMATRIX( Q, c, 7 ) = /* error */
- YMATRIX( Q, c, 6 ) / YMATRIX( Q, c, 4 ) ;
+ /* Compute R = AX - B */
+ gsl_vector_memcpy(R, AX);
+ gsl_vector_sub( R, B );
+
+ Q = gsl_matrix_alloc( count, 7 );
+ for( c = 0; c < count; c++ ){
+ gsl_matrix_set(Q, c, 0, gsl_matrix_get(A, c, 1)) ; /* x */
+ gsl_matrix_set(Q, c, 1, gsl_matrix_get(A, c, 3)) ; /* y */
+ gsl_matrix_set(Q, c, 2, gsl_matrix_get(A, c, 5)) ; /* pins */
+ gsl_matrix_set(Q, c, 3, gsl_vector_get(B, c)) ; /* B */
+ gsl_matrix_set(Q, c, 4, gsl_matrix_get(AX, c, 0)) ; /* AX */
+ gsl_matrix_set(Q, c, 5, gsl_matrix_get(R, c, 0)) ; /* B - AX */
+ gsl_matrix_set(Q, c, 6, gsl_matrix_get(Q, c, 5) /
+ gsl_matrix_get(Q, c, 3)) ; /* error */
}
fprintf( stderr, "\n x, y, B, pins, AB, B - AX, error:\n");
- Ymatrix_disp( Q ) ;
+ gsl_matrix_disp( Q, count, 7 ) ;
+ gsl_matrix_free(Q);
+ gsl_vector_free(R);
+ gsl_vector_free(AX);
) ;
/* now done free the matrices */
- Ymatrix_free( A ) ;
- Ymatrix_free( B ) ;
- Ymatrix_free( Xret ) ;
+ gsl_matrix_free( A );
+ gsl_vector_free( B );
+ gsl_vector_free( Xret );
} /* end adapt_wire_estimator */
diff -Nru graywolf-0.1.4/src/Ylib/CMakeLists.txt graywolf-0.1.4+20170306gitecee764/src/Ylib/CMakeLists.txt
--- graywolf-0.1.4/src/Ylib/CMakeLists.txt 2016-10-06 22:10:19.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/src/Ylib/CMakeLists.txt 2017-03-06 20:47:22.000000000 +0100
@@ -3,7 +3,7 @@
#add_executable(mincut main.c output.c readcells.c ${CMAKE_SOURCE_DIR}/src/date/date.c)
-add_library(ycadgraywolf SHARED assign.c buster.c cleanup.c colors.c deck.c dialog.c draw.c dset.c edcolors.c file.c getftime.c graph.c grid.c hash.c heap.c list.c log.c matrix.c menus.c message.c mst.c mytime.c okmalloc.c path.c plot.c program.c project.c queue.c quicksort.c radixsort.c rand.c rbtree.c relpath.c set.c stat.c stats.c string.c svd.c system.c time.c timer.c trans.c wgraphics.c ydebug.c yreadpar.c )
+add_library(ycadgraywolf SHARED assign.c buster.c cleanup.c colors.c deck.c dialog.c draw.c dset.c edcolors.c file.c getftime.c graph.c grid.c hash.c heap.c list.c log.c menus.c message.c mst.c mytime.c okmalloc.c path.c plot.c program.c project.c queue.c quicksort.c radixsort.c rand.c rbtree.c relpath.c set.c stat.c stats.c string.c system.c time.c timer.c trans.c wgraphics.c ydebug.c yreadpar.c )
target_link_libraries(ycadgraywolf X11)
target_link_libraries(ycadgraywolf m)
diff -Nru graywolf-0.1.4/src/Ylib/matrix.c graywolf-0.1.4+20170306gitecee764/src/Ylib/matrix.c
--- graywolf-0.1.4/src/Ylib/matrix.c 2016-10-06 22:10:19.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/src/Ylib/matrix.c 1970-01-01 01:00:00.000000000 +0100
@@ -1,451 +0,0 @@
-/*
- * Copyright (C) 1991 Yale University
- *
- * This work is distributed in the hope that it will be useful; 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 any later version, on the following conditions:
- *
- * (a) YALE MAKES NO, AND EXPRESSLY DISCLAIMS
- * ALL, REPRESENTATIONS OR WARRANTIES THAT THE MANUFACTURE, USE, PRACTICE,
- * SALE OR
- * OTHER DISPOSAL OF THE SOFTWARE DOES NOT OR WILL NOT INFRINGE UPON ANY
- * PATENT OR
- * OTHER RIGHTS NOT VESTED IN YALE.
- *
- * (b) YALE MAKES NO, AND EXPRESSLY DISCLAIMS ALL, REPRESENTATIONS AND
- * WARRANTIES
- * WHATSOEVER WITH RESPECT TO THE SOFTWARE, EITHER EXPRESS OR IMPLIED,
- * INCLUDING,
- * BUT NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
- * PARTICULAR
- * PURPOSE.
- *
- * (c) LICENSEE SHALL MAKE NO STATEMENTS, REPRESENTATION OR WARRANTIES
- * WHATSOEVER TO
- * ANY THIRD PARTIES THAT ARE INCONSISTENT WITH THE DISCLAIMERS BY YALE IN
- * ARTICLE
- * (a) AND (b) above.
- *
- * (d) IN NO EVENT SHALL YALE, OR ITS TRUSTEES, DIRECTORS, OFFICERS,
- * EMPLOYEES AND
- * AFFILIATES BE LIABLE FOR DAMAGES OF ANY KIND, INCLUDING ECONOMIC DAMAGE OR
- * INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER YALE SHALL BE
- * ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE
- * POSSIBILITY OF THE FOREGOING.
- *
- */
-
-/* -----------------------------------------------------------------
-FILE: matrix.c
-DESCRIPTION:a curvefitting program.
-CONTENTS:
-DATE: Tue Jan 15 01:35:11 EST 1991 - original coding.
-REVISIONS: Sun Nov 3 12:48:39 EST 1991 - made matrix memory
- allocation more efficient using YVECTOR routines.
- 12/09/91 - Add prototype macro for non ansi compiler -R.A.Weier
------------------------------------------------------------------ */
-#ifndef lint
-static char SccsId[] = "@(#) matrix.c version 1.3 12/9/91" ;
-#endif
-
-#include <yalecad/base.h>
-#include <yalecad/message.h>
-#include <yalecad/debug.h>
-#include <yalecad/linalg.h>
-
-static DOUBLE find_det( P3(DOUBLE **a, INT rows, INT columns ) );
-static sign_cof( P2(INT row, INT column) ) ;
-
-/* *****************************************************************
- Allocate the space for the matrix.
- **************************************************************** */
-YMPTR Ymatrix_create( rows, columns )
-INT rows, columns ;
-{
- YMPTR mptr ; /* pointer to matrix record */
- DOUBLE **m ; /* pointer to matrix space */
- INT i ;
-
- mptr = YMALLOC( 1, YMBOX ) ;
- /* now allocate memory for rows */
- m = mptr->m = YVECTOR_MALLOC( 1, rows, DOUBLE * ) ;
- /* now allocate memory for columns */
- for( i = 1;i<= rows; i++ ){
- m[i] = YVECTOR_CALLOC( 1, columns, DOUBLE ) ;
- }
- mptr->rows = rows ;
- mptr->columns = columns ;
- return( mptr ) ;
-} /* end Ymatrix_create */
-
-/* *****************************************************************
- Free the space of a matrix.
- **************************************************************** */
-YMPTR Ymatrix_free( mptr )
-YMPTR mptr ; /* pointer to matrix record */
-{
- DOUBLE **m ; /* pointer to matrix space */
- INT i ;
- INT rows ;
-
- rows = mptr->rows ;
- m = mptr->m ;
- /* free memory */
- for( i = 1; i <= rows; i++ ){
- YVECTOR_FREE( m[i], 1 ) ;
- }
- YVECTOR_FREE( m, 1 ) ;
- YFREE( mptr ) ;
-} /* end Ymatrix_free */
-
-YMPTR Ymatrix_transpose( mptr )
-YMPTR mptr ;
-{
- INT i, j ;
- DOUBLE **m ;
- DOUBLE **b ;
- YMPTR buf ;
-
- buf = Ymatrix_create( mptr->columns, mptr->rows ) ;
-
- m = mptr->m ;
- b = buf->m ;
-
- for( i=1; i <= mptr->rows; i++ ){
- for( j=1; j <= mptr->columns; j++ ){
- b[j][i] = m[i][j] ;
- }
- }
- return( buf ) ;
-
-} /* end Ymatrix_transpose */
-
-YMPTR Ymatrix_mult( aptr, bptr )
-YMPTR aptr, bptr ;
-{
- INT i, j, k, n ;
- DOUBLE **a ;
- DOUBLE **b ;
- DOUBLE **c ;
- DOUBLE result ;
- YMPTR buf ;
-
- /* first perform error checking */
- if( aptr->columns != bptr->rows ){
- sprintf( YmsgG,
- "Matrices cannot be multiplied together A (%dx%d) B(%dx%d)\n",
- aptr->rows, aptr->columns, bptr->rows, bptr->columns ) ;
- M( ERRMSG, "Ymatrix_mult", YmsgG ) ;
- }
- buf = Ymatrix_create( aptr->rows, bptr->columns ) ;
-
- a = aptr->m ;
- b = bptr->m ;
- c = buf->m ;
- n = aptr->columns ;
-
- for( i=1; i <= aptr->rows; i++ ){
- for( j=1; j <= bptr->columns; j++ ){
- result = 0.0 ;
- for( k=1; k<= n; k++ ){
- result += a[i][k] * b[k][j] ;
- }
- c[i][j] = result ;
- }
- }
- return( buf ) ;
-
-} /* end Ymatrix_mult */
-
-YMPTR Ymatrix_sub( aptr, bptr )
-YMPTR aptr, bptr ;
-{
- INT i, j ;
- DOUBLE **a, **b, **c ;
- DOUBLE *fast_a, *fast_b, *fast_c ;
- DOUBLE result ;
- YMPTR buf ;
-
- /* first perform error checking */
- if( aptr->rows != bptr->rows ||
- aptr->columns != bptr->columns ){
- sprintf( YmsgG,
- "Matrices cannot be subtracted A (%dx%d) B(%dx%d)\n",
- aptr->rows, aptr->columns, bptr->rows, bptr->columns ) ;
- M( ERRMSG, "Ymatrix_sub", YmsgG ) ;
- }
- buf = Ymatrix_create( aptr->rows, aptr->columns ) ;
-
- a = aptr->m ;
- b = bptr->m ;
- c = buf->m ;
- for( i=1; i <= aptr->rows; i++ ){
- fast_a = a[i] ;
- fast_b = b[i] ;
- fast_c = c[i] ;
- for( j=1; j <= aptr->columns; j++ ){
- fast_c[j] = fast_a[j] - fast_b[j] ;
- }
- }
- return( buf ) ;
-
-} /* end Ymatrix_sub */
-
-Ymatrix_disp( mptr )
-YMPTR mptr ;
-{
- DOUBLE **m ;
- INT i, j ;
- m = mptr->m ;
-
- for( i=1; i <= mptr->rows; i++ ){
- for( j=1; j <= mptr->columns; j++ ){
- fprintf( stderr, "% 4.4le ", m[i][j] ) ;
- }
- fprintf( stderr, "\n" ) ;
- }
- fprintf( stderr, "\n" ) ;
-} /* end Ymatrix_disp */
-
-YMPTR Ymatrix_eye( size )
-INT size ;
-{
- INT i, j ;
- DOUBLE **m ;
- DOUBLE *fast ;
- YMPTR buf ;
-
- buf = Ymatrix_create( size, size ) ;
- m = buf->m ;
- for( i = 1; i <= size; i++ ){
- fast = m[i] ;
- for( j = 1; j <= size ; j++ ){
- if( i == j ){
- /* one is exact in fp. */
- fast[j] = 1.0 ;
- } else {
- fast[j] = YZERO ;
- }
- }
- }
- return( buf ) ;
-} /* end Ymatrix_eye */
-
-Ymatrix_zero( matrix )
-YMPTR matrix ;
-{
- INT r, c, i, j ;
- DOUBLE *fast ;
- DOUBLE **m ;
-
- m = matrix->m ;
- r = matrix->rows ;
- c = matrix->columns ;
- for( i = 1; i <= r; i++ ){
- fast = m[i] ;
- for( j = 1; j <= c ; j++ ){
- fast[j] = YZERO ;
- }
- }
-} /* end Ymatrix_zero */
-
-
-YMPTR Ymatrix_copy( input )
-YMPTR input ;
-{
- INT i, j ;
- DOUBLE **in_mat, **copy_mat ;
- DOUBLE *in_fast, *copy_fast ;
- YMPTR copy ;
-
- copy = Ymatrix_create( input->rows, input->columns ) ;
- in_mat = input->m ;
- copy_mat = copy->m ;
- for( i = 1; i <= input->rows; i++ ){
- in_fast = in_mat[i] ;
- copy_fast = copy_mat[i] ;
- for( j = 1; j <= input->columns ; j++ ){
- copy_fast[j] = in_fast[j] ;
- }
- }
- return( copy ) ;
-} /* end Ymatrix_copy */
-
-YMPTR Ymatrix_linv( aptr )
-YMPTR aptr ;
-{
- INT i, j, k, n ;
- DOUBLE **a ;
- DOUBLE **b ;
- DOUBLE **c ;
- DOUBLE det, recip_det ;
- YMPTR cof ;
- YMPTR buf ;
-
- buf = Ymatrix_create( aptr->rows, aptr->columns ) ;
- cof = Ymatrix_cofactors( aptr ) ;
-
- det = 0.0 ;
- a = aptr->m ;
- c = cof->m ;
- /* find the determinant using first column */
- for( i=1; i <= aptr->rows; i++ ){
- det += a[i][1] * c[i][1] ;
- }
-
- ASSERT( det == find_det( aptr->m, aptr->rows, aptr->columns ),
- "Yinv", "Problem with determinant" ) ;
-
- if( det != 0.0 ){
- recip_det = 1.0 / det ;
- buf = Ymatrix_transpose( cof ) ;
- b = buf->m ;
- for( i=1; i <= buf->rows; i++ ){
- for( j=1; j <= buf->columns; j++ ){
- b[i][j] *= recip_det ;
- }
- }
- }
- return( buf ) ;
-
-} /* end Ymatrix_linv */
-
-YMPTR Ymatrix_cofactors( aptr )
-YMPTR aptr ;
-{
-
- INT i, j, k, l ;
- INT r, c ;
- INT rows, columns ;
- DOUBLE cofactor ;
- DOUBLE **a, **b, **m ;
- YMPTR buf ;
-
- buf = Ymatrix_create( aptr->rows, aptr->columns ) ;
- rows = aptr->rows ;
- columns = aptr->columns ;
- a = aptr->m ; /* a matrix - want cofactors of this matrix */
- b = buf->m ; /* buf matrix */
-
- /* allocate space for a row - 1 x column - 1 clone of a matrix */
- m = YVECTOR_MALLOC( 1, rows - 1, DOUBLE * ) ;
- /* now allocate memory for columns */
- for( i = 1;i< rows; i++ ){
- m[i] = YVECTOR_CALLOC( 1, columns-1, DOUBLE ) ;
- }
-
- /* process all cofactors */
- for( i=1; i <= rows; i++ ){
- for( j=1; j <= columns; j++ ){
-
- /* load cofactor array */
- /* strike row and column */
- r = 0 ;
- for( k = 1; k <= rows; k++ ){
- if( k == i ){
- continue ;
- }
- r++ ;
- c = 0 ;
- for( l = 1; l <= columns; l++ ){
- if( l == j ){
- continue ;
- }
- m[r][++c] = a[k][l] ;
- }
- } /* end load of cofactor array */
-
- /* process one cofactor */
- cofactor = sign_cof(i,j) * find_det( m, rows-1, columns-1 ) ;
- /* store result in buf */
- b[i][j] = cofactor ;
-
- }
- }
- /* now time to free memory */
- for( i = 1;i< columns; i++ ){
- YVECTOR_FREE( m[i], 1 ) ;
- }
- YVECTOR_FREE( m, 1 ) ;
- return( buf ) ;
-} /* end Ymatrix_cofactors */
-
-/* find a determinant - works on the actual memory arrays - not user */
-/* records */
-static DOUBLE find_det( a, rows, columns )
-DOUBLE **a ;
-INT rows, columns ;
-{
- DOUBLE result ;
- DOUBLE cofactor ;
- DOUBLE **m ;
- INT i, j, k ;
- INT r, c ;
-
- if( rows == 1 && columns == 1 ){
- result = a[1][1] ;
-
- } else if( rows == 2 && columns == 2 ){
- result = a[1][1] * a[2][2] - a[2][1] * a[1][2] ;
-
- } else {
-
- /* allocate space for an row-1 x columns-1 matrix */
- m = YVECTOR_MALLOC( 1, rows-1, DOUBLE * ) ;
- /* now allocate memory for columns */
- for( i = 1;i< rows; i++ ){
- m[i] = YVECTOR_MALLOC( 1, columns-1, DOUBLE ) ;
- }
-
- /* do Laplace expansion along top row - all columns */
- result = 0.0 ;
- for( k=1; k <= columns; k++ ){
-
- /* load array by striking row and column */
- /* struck row is always one */
- r = 0 ;
- for( i = 2; i <= rows; i++ ){
- r++ ;
- c = 0 ;
- for( j = 1; j <= columns; j++ ){
- if( j == k ){
- continue ;
- }
- m[r][++c] = a[i][j] ;
- }
- } /* end load of cofactor array */
-
- ASSERT( r == rows-1, "Yfind_det","problem with rows" ) ;
- ASSERT( c == columns-1, "Yfind_det","problem with rows" ) ;
-
- /* recursively call find_det to get result */
- cofactor = sign_cof(1,k) * find_det( m, rows-1, columns-1 ) ;
- result += a[1][k] * cofactor ;
- }
-
- /* now time to free memory */
- for( i = 1;i< columns; i++ ){
- YVECTOR_FREE( m[i], 1 ) ;
- }
- YVECTOR_FREE( m, 1 ) ;
- }
-
- return( result ) ;
-
-} /* end find_det */
-
-static sign_cof( row, column )
-INT row, column ;
-{
- INT sum ;
-
- /* implements (-1) ** (row + column) */
- sum = row + column ;
- if( sum % 2 ){ /* odd power */
- return( -1.0 ) ;
- } else { /* even power */
- return( 1.0 ) ;
- }
-} /* end sign_cof */
diff -Nru graywolf-0.1.4/src/Ylib/svd.c graywolf-0.1.4+20170306gitecee764/src/Ylib/svd.c
--- graywolf-0.1.4/src/Ylib/svd.c 2016-10-06 22:10:19.000000000 +0200
+++ graywolf-0.1.4+20170306gitecee764/src/Ylib/svd.c 1970-01-01 01:00:00.000000000 +0100
@@ -1,522 +0,0 @@
-/*
- * Copyright (C) 1991 Yale University
- *
- * This work is distributed in the hope that it will be useful; 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 any later version, on the following conditions:
- *
- * (a) YALE MAKES NO, AND EXPRESSLY DISCLAIMS
- * ALL, REPRESENTATIONS OR WARRANTIES THAT THE MANUFACTURE, USE, PRACTICE,
- * SALE OR
- * OTHER DISPOSAL OF THE SOFTWARE DOES NOT OR WILL NOT INFRINGE UPON ANY
- * PATENT OR
- * OTHER RIGHTS NOT VESTED IN YALE.
- *
- * (b) YALE MAKES NO, AND EXPRESSLY DISCLAIMS ALL, REPRESENTATIONS AND
- * WARRANTIES
- * WHATSOEVER WITH RESPECT TO THE SOFTWARE, EITHER EXPRESS OR IMPLIED,
- * INCLUDING,
- * BUT NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
- * PARTICULAR
- * PURPOSE.
- *
- * (c) LICENSEE SHALL MAKE NO STATEMENTS, REPRESENTATION OR WARRANTIES
- * WHATSOEVER TO
- * ANY THIRD PARTIES THAT ARE INCONSISTENT WITH THE DISCLAIMERS BY YALE IN
- * ARTICLE
- * (a) AND (b) above.
- *
- * (d) IN NO EVENT SHALL YALE, OR ITS TRUSTEES, DIRECTORS, OFFICERS,
- * EMPLOYEES AND
- * AFFILIATES BE LIABLE FOR DAMAGES OF ANY KIND, INCLUDING ECONOMIC DAMAGE OR
- * INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER YALE SHALL BE
- * ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE
- * POSSIBILITY OF THE FOREGOING.
- *
- */
-
-/* -----------------------------------------------------------------
-FILE: svd.c
-DESCRIPTION:Singular value decomposition routines.
-CONTENTS:
-DATE: Tue Jan 15 00:52:31 EST 1991 - original coding.
-REVISIONS: Tue Jan 15 22:26:57 PST 1991 - added TRUE return condition.
- Thu Jan 17 00:59:42 PST 1991 - now notify user of
- singular value when in debug mode.
- Thu Jan 24 20:17:23 PST 1991 - changed vector routines.
- Tue Mar 12 16:56:22 CST 1991 - got rid of unnecessary
- assignment.
------------------------------------------------------------------ */
-#ifndef lint
-static char SccsId[] = "@(#) svd.c version 1.7 12/15/91" ;
-#endif
-
-#include <yalecad/base.h>
-#include <yalecad/message.h>
-#include <yalecad/debug.h>
-#include <yalecad/linalg.h>
-
-/* static definitions */
-static DOUBLE atS,btS,ctS;
-static DOUBLE maxarg1S,maxarg2S;
-
-#undef MAX
-#define PYTHAG(a,b) ((atS=fabs(a)) > (btS=fabs(b)) ? \
-(ctS=btS/atS,atS*sqrt(1.0+ctS*ctS)) : (btS ? (ctS=atS/btS,btS*sqrt(1.0+ctS*ctS)): 0.0))
-#define MAX(a,b) (maxarg1S=(a),maxarg2S=(b),(maxarg1S) > (maxarg2S) ?\
- (maxarg1S) : (maxarg2S))
-#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
-#define TOLERANCE 1e-6
-
-/* ---------------------------------------------------------------
- Procedure Ysvd_solve - solve a set of linear equations using SVD.
- Input: A - a m x n matrix.
- Input: B - a m x k matrix.
- Output: X - a n x k matrix.
------------------------------------------------------------------*/
-BOOL Ysvd_solve( A, B, Xret )
-YMPTR A, B, *Xret ;
-{
- INT i, j, jj, k, m, n, t ;
- YMPTR U, W, V ;
- DOUBLE wmax, threshold, sum ;
- DOUBLE **b, **u, **v, **w, **x ;
- DOUBLE *tmp ;
-
- m = A->rows ;
- n = A->columns ;
- k = B->columns ;
-
- /* next make sure B conforms to A */
- if( m != B->rows ){
- M( ERRMSG, "Ysvd_solve", "The matrices A and B do not conform\n" ) ;
- sprintf( YmsgG, "A = %d x %d B = %d x %d\n\n", m, n, B->rows, k ) ;
- M( ERRMSG, NULL, YmsgG ) ;
- return(FALSE) ;
- }
-
- /* first factor A using SVD */
- if(!(Ysvd_decompose( A, &U, &W, &V ))){
- return( FALSE ) ;
- }
- D( "Ysvd_solve/decompose",
- YMPTR Vt ;
- YMPTR WVt ;
- YMPTR UWVt ;
-
- fprintf( stderr, "U:\n" ) ;
- Ymatrix_disp( U ) ;
- fprintf( stderr, "W:\n" ) ;
- Ymatrix_disp( W ) ;
- fprintf( stderr, "V:\n" ) ;
- Ymatrix_disp( V ) ;
-
- Vt = Ymatrix_transpose( V ) ;
- fprintf( stderr, "Vt:\n" ) ;
- Ymatrix_disp( Vt ) ;
- WVt = Ymatrix_mult( W, Vt ) ;
- UWVt = Ymatrix_mult( U, WVt ) ;
- fprintf( stderr, "U * W * Vt:\n" ) ;
- Ymatrix_disp( UWVt ) ;
- ) ;
-
- /* now throw away those things which make the answer unstable */
- /* first find the maximum w */
- wmax=0.0;
- w = W->m ;
- for (j=1;j<=n;j++){
- if(w[j][j] > wmax){
- wmax=w[j][j] ;
- }
- }
- /* next create a threshold */
- threshold = TOLERANCE * wmax ;
- /* now zero below threshold */
- for (j=1;j<=n;j++){
- if( w[j][j] < threshold ){
- D( "Ysvd_solve/singular",
- fprintf( stderr, "Singular value for w[%d]:%4.6le\n",
- j, w[j][j] ) ;
- ) ;
- w[j][j] = 0.0;
- }
- }
- D( "Ysvd_solve/a_weight",
- fprintf( stderr, "W:\n" ) ;
- Ymatrix_disp( W ) ;
- ) ;
-
- D( "Ysvd_solve/backsub",
- YMPTR Ut ;
- YMPTR Utb ;
- YMPTR WUtb ;
- YMPTR Xver ;
- YMPTR Wprime ;
- DOUBLE **wp ;
-
- Ut = Ymatrix_transpose( U ) ;
- Utb = Ymatrix_mult( Ut, B ) ;
- Wprime = Ymatrix_copy( W ) ;
- wp = Wprime->m ;
- for( i = 1; i <= n; i++ ){
- if( wp[i][i] ){
- wp[i][i] = 1 / wp[i][i] ;
- }
- }
- WUtb = Ymatrix_mult( Wprime, Utb ) ;
- Xver = Ymatrix_mult( V, WUtb ) ;
- fprintf( stderr, "V * 1 / w * Ut * B:\n" ) ;
- Ymatrix_disp( Xver ) ;
- ) ;
-
- /* now we are ready for back substitution */
- /* x = V * diag( 1 / w ) * Ut * b */
- *Xret = Ymatrix_create( n, k ) ;
- x = (*Xret)->m ;
- u = U->m ;
- v = V->m ;
- b = B->m ;
- tmp = YVECTOR_MALLOC(1,n,DOUBLE) ;
- for( t = 1; t <= k ; t++ ){
- /* first calculate Ut * b */
- for (j=1;j<=n;j++) {
- sum=0.0;
- if (w[j][j]) {
- for (i=1;i<=m;i++){
- sum += u[i][j]*b[i][t];
- }
- /* this is the divide by 1 / w */
- sum /= w[j][j];
- }
- tmp[j]=sum;
- }
- /* now multiply by V */
- for (j=1;j<=n;j++) {
- sum=0.0;
- for (jj=1;jj<=n;jj++){
- sum += v[j][jj]*tmp[jj];
- }
- x[j][t]=sum;
- }
- } /* end looping on columns */
-
- /* free memory */
- YVECTOR_FREE(tmp,1) ;
- Ymatrix_free( U ) ;
- Ymatrix_free( V ) ;
- Ymatrix_free( W ) ;
- return( TRUE ) ;
-
-} /* end Ysvd_solve */
-
-/* -------------------------------------------------------------
- Procedure Ysvd_deccompose - decompose A into U W Vt.
- Input: A an m x n matrix.
- Output: U - a m x n column orthogonal matrix.
- Output: W - a n x 1 column vector.
- Output: V - a n x n column orthogonal matrix.
----------------------------------------------------------------- */
-BOOL Ysvd_decompose( A, Uret, Wret, Vret )
-YMPTR A, *Uret, *Wret, *Vret ;
-{
- INT m, n ;
- INT flag,i,its,j,jj,k,l,nm;
- DOUBLE c,f,h,s,x,y,z;
- DOUBLE anorm=0.0,g=0.0,scale=0.0;
- DOUBLE *rv1;
- DOUBLE **a,*w,**v, **wfast ;
-
- m = A->rows ;
- n = A->columns ;
-
- if (m < n){
- sprintf( YmsgG, "m < n (A = %d x %d)\n", m, n ) ;
- M( ERRMSG, "Ysvd_decompose", YmsgG ) ;
- M( ERRMSG, NULL, "You must augment A with extra zero rows");
- return( FALSE ) ;
- }
-
- /* create the working space */
- *Uret = Ymatrix_copy( A ) ;
- a = (*Uret)->m ;
- *Wret = Ymatrix_eye( n ) ;
- *Vret = Ymatrix_create( n, n ) ;
- v = (*Vret)->m ;
- w = YVECTOR_MALLOC(1,n,DOUBLE);
- rv1 = YVECTOR_MALLOC(1,n,DOUBLE);
-
- /* Householder reduction to bidiagonal form */
- for (i=1;i<=n;i++) {
- l=i+1;
- rv1[i]=scale*g;
- g=s=scale=0.0;
- if (i <= m) {
- for (k=i;k<=m;k++) scale += fabs(a[k][i]);
- if (scale) {
- for (k=i;k<=m;k++) {
- a[k][i] /= scale;
- s += a[k][i]*a[k][i];
- }
- f=a[i][i];
- g = -SIGN(sqrt(s),f);
- h=f*g-s;
- a[i][i]=f-g;
- if (i != n) {
- for (j=l;j<=n;j++) {
- for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
- f=s/h;
- for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
- }
- }
- for (k=i;k<=m;k++) a[k][i] *= scale;
- }
- }
- w[i]=scale*g;
- g=s=scale=0.0;
- if (i <= m && i != n) {
- for (k=l;k<=n;k++) scale += fabs(a[i][k]);
- if (scale) {
- for (k=l;k<=n;k++) {
- a[i][k] /= scale;
- s += a[i][k]*a[i][k];
- }
- f=a[i][l];
- g = -SIGN(sqrt(s),f);
- h=f*g-s;
- a[i][l]=f-g;
- for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
- if (i != m) {
- for (j=l;j<=m;j++) {
- for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
- for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
- }
- }
- for (k=l;k<=n;k++) a[i][k] *= scale;
- }
- }
- anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i])));
- } /* end of Householder reduction */
-
- /* Accummulation of right-hand transformations */
- for (i=n;i>=1;i--) {
- if (i < n) {
- if (g) {
- for (j=l;j<=n;j++)
- v[j][i]=(a[i][j]/a[i][l])/g;
- for (j=l;j<=n;j++) {
- for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
- for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
- }
- }
- for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
- }
- v[i][i]=1.0;
- g=rv1[i];
- l=i;
- } /* end accummulation of right-hand transformation */
-
- /* Accummulation of left-hand transformations */
- for (i=n;i>=1;i--) {
- l=i+1;
- g=w[i];
- if (i < n) for (j=l;j<=n;j++) a[i][j]=0.0;
- if (g) {
- g=1.0/g;
- if (i != n) {
- for (j=l;j<=n;j++) {
- for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
- f=(s/a[i][i])*g;
- for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
- }
- }
- for (j=i;j<=m;j++) a[j][i] *= g;
- } else {
- for (j=i;j<=m;j++) a[j][i]=0.0;
- }
- ++a[i][i];
- } /* end accummulation of left-hand transformation */
-
- /* Diagonalization of the bidiagonal form */
- /* loop over singular values */
- for (k=n;k>=1;k--) {
- /* loop over allowed iterations */
- for (its=1;its<=30;its++) {
- flag=1;
- for (l=k;l>=1;l--) {
- nm=l-1;
- if (fabs(rv1[l])+anorm == anorm) {
- flag=0;
- break;
- }
- if (fabs(w[nm])+anorm == anorm) break;
- }
- if (flag) {
- s=1.0;
- for (i=l;i<=k;i++) {
- f=s*rv1[i];
- if (fabs(f)+anorm != anorm) {
- g=w[i];
- h=PYTHAG(f,g);
- w[i]=h;
- h=1.0/h;
- c=g*h;
- s=(-f*h);
- for (j=1;j<=m;j++) {
- y=a[j][nm];
- z=a[j][i];
- a[j][nm]=y*c+z*s;
- a[j][i]=z*c-y*s;
- }
- }
- }
- }
- /* Convergence */
- z=w[k];
- if (l == k) {
- if (z < 0.0) {
- w[k] = -z;
- for (j=1;j<=n;j++) v[j][k]=(-v[j][k]);
- }
- break;
- }
- if (its == 30){
- M(ERRMSG, "Ysvd_decompose", "No convergence in 30 SVDCMP iterations");
- return( FALSE ) ;
- }
- x=w[l];
- nm=k-1;
- y=w[nm];
- g=rv1[nm];
- h=rv1[k];
- f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
- g=PYTHAG(f,1.0);
- f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
-
- /* The QR transformation */
- c=s=1.0;
- for (j=l;j<=nm;j++) {
- i=j+1;
- g=rv1[i];
- y=w[i];
- h=s*g;
- g=c*g;
- z=PYTHAG(f,h);
- rv1[j]=z;
- c=f/z;
- s=h/z;
- f=x*c+g*s;
- g=g*c-x*s;
- h=y*s;
- y=y*c;
- for (jj=1;jj<=n;jj++) {
- x=v[jj][j];
- z=v[jj][i];
- v[jj][j]=x*c+z*s;
- v[jj][i]=z*c-x*s;
- }
- z=PYTHAG(f,h);
- w[j]=z;
- if (z) {
- z=1.0/z;
- c=f*z;
- s=h*z;
- }
- f=(c*g)+(s*y);
- x=(c*y)-(s*g);
- for (jj=1;jj<=m;jj++) {
- y=a[jj][j];
- z=a[jj][i];
- a[jj][j]=y*c+z*s;
- a[jj][i]=z*c-y*s;
- }
- }
- rv1[l]=0.0;
- rv1[k]=f;
- w[k]=x;
- }
- }
-
- /* now transfer w to Wret */
- wfast = (*Wret)->m ;
- for( i = 1; i <= n; i++ ){
- wfast[i][i] = w[i] ;
- }
- YVECTOR_FREE(rv1,1) ;
- YVECTOR_FREE(w,1) ;
- return( TRUE ) ;
-
-} /* end Ysvd_decompose */
-
-#undef SIGN
-#undef MAX
-#undef PYTHAG
-
-#ifdef TEST
-/* *************************** TEST SOLVER ****************** */
-#include <yalecad/cleanup.h>
-main(argc,argv)
-int argc;
-char *argv[];
-{
-
- YMPTR X, y, Xt, XtX, B, Xret ;
-
- /* start up cleanup handler */
- YINITCLEANUP( argv[0], NULL, MAYBEDUMP ) ;
- YsetDebug( TRUE ) ;
-
- X = Ymatrix_create( 5, 3 ) ;
- y = Ymatrix_create( 5, 1 ) ;
-
- /* load matrix */
- YMATRIX( X, 1, 1 ) = 1.0 ;
- YMATRIX( X, 1, 2 ) = 0.0 ;
- YMATRIX( X, 1, 3 ) = 0.0 ;
- YMATRIX( X, 2, 1 ) = 1.0 ;
- YMATRIX( X, 2, 2 ) = 67.0 ;
- YMATRIX( X, 2, 3 ) = 448.0 ;
- YMATRIX( X, 3, 1 ) = 1.0 ;
- YMATRIX( X, 3, 2 ) = 157.0 ;
- YMATRIX( X, 3, 3 ) = 24649.0 ;
- YMATRIX( X, 4, 1 ) = 1.0 ;
- YMATRIX( X, 4, 2 ) = 67.0 ;
- YMATRIX( X, 4, 3 ) = 448.0 ;
- YMATRIX( X, 5, 1 ) = 1.0 ;
- YMATRIX( X, 5, 2 ) = 0.0 ;
- YMATRIX( X, 5, 3 ) = 0.0 ;
- M( MSG, NULL, "X:\n" ) ;
- Ymatrix_disp( X ) ;
-
- YMATRIX( y, 1, 1 ) = 0.7 ;
- YMATRIX( y, 2, 1 ) = 0.7 ;
- YMATRIX( y, 3, 1 ) = 0.15 ;
- YMATRIX( y, 4, 1 ) = 0.7 ;
- YMATRIX( y, 5, 1 ) = 0.7 ;
- M( MSG, NULL, "y:\n" ) ;
- Ymatrix_disp( y ) ;
-
- M( MSG, NULL, "\nXt:\n" ) ;
- Xt = Ymatrix_transpose( X ) ;
- Ymatrix_disp( Xt ) ;
-
- M( MSG, NULL, "\nXtX:\n" ) ;
- XtX = Ymatrix_mult( Xt, X ) ;
- Ymatrix_disp( XtX ) ;
-
- /* now test eq solver */
- B = Ymatrix_mult( Xt, y ) ;
- M( MSG, NULL, "\nB originally\n");
- Ymatrix_disp( B ) ;
-
- if( Ysvd_solve( XtX, B, &Xret )){
- M( MSG, NULL, "\nSvd Solver\n");
- Ymatrix_disp( Xret ) ;
- }
-
- YexitPgm(0) ;
-
-} /* end main program */
-
-#endif /* TEST */
--- End Message ---