diff --git a/CMakeLists.txt b/CMakeLists.txt index 10e2b22..4bc06e5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 --git a/debian/changelog b/debian/changelog index a375b13..3126bc9 100644 --- a/debian/changelog +++ b/debian/changelog @@ -1,3 +1,29 @@ +graywolf (0.1.4+20170307gite1bf319-2) unstable; urgency=medium + + * debian/copyright: + - Added Staf Verhaegen (should have been added before uploading 0.1.4-1) + (Closes: #857088) + - Added more copyright years for Tim Edwards and Ruben Undheim + + -- Ruben Undheim Wed, 08 Mar 2017 09:25:41 +0100 + +graywolf (0.1.4+20170307gite1bf319-1) unstable; urgency=medium + + * New upstream GIT HEAD with bumped SOVERSION and major version + * New binary package libycadgraywolf1 + + -- Ruben Undheim Tue, 07 Mar 2017 20:40:23 +0000 + +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 Mon, 06 Mar 2017 20:26:57 +0000 + graywolf (0.1.4-1) unstable; urgency=low * New upstream release diff --git a/debian/control b/debian/control index e73b0e1..94a2cb7 100644 --- a/debian/control +++ b/debian/control @@ -5,13 +5,15 @@ Section: electronics 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 Homepage: https://github.com/rubund/graywolf -Package: libycadgraywolf0 +Package: libycadgraywolf1 Architecture: any Section: libs Depends: ${shlibs:Depends}, @@ -27,7 +29,7 @@ Description: Library for common functions used by graywolf Package: libycadgraywolf-dev Architecture: any Section: libdevel -Depends: libycadgraywolf0 (= ${binary:Version}), +Depends: libycadgraywolf1 (= ${binary:Version}), ${misc:Depends} Multi-Arch: same Description: Library for graywolf - development files @@ -36,14 +38,14 @@ Description: Library for graywolf - development files . This is the development files for the library. -Package: libycadgraywolf0-dbg +Package: libycadgraywolf1-dbg Architecture: any Section: debug Priority: extra -Depends: libycadgraywolf0 (= ${binary:Version}), +Depends: libycadgraywolf1 (= ${binary:Version}), ${misc:Depends} Multi-Arch: same -Description: debug symbols for libycadgraywolf0 +Description: debug symbols for libycadgraywolf1 This library contains functions shared by the graywolf binaries found in the /usr/lib/graywolf/bin/ directory. . @@ -51,7 +53,7 @@ Description: debug symbols for libycadgraywolf0 Package: graywolf Architecture: any -Depends: libycadgraywolf0 (= ${binary:Version}), +Depends: libycadgraywolf1 (= ${binary:Version}), ${shlibs:Depends}, ${misc:Depends} Recommends: qflow diff --git a/debian/copyright b/debian/copyright index 5ccd989..eb6e6bb 100644 --- a/debian/copyright +++ b/debian/copyright @@ -4,8 +4,9 @@ Source: https://github.com/rubund/graywolf Files: * Copyright: 1987-1992 Yale University - 2006-2015 Tim Edwards - 2014 Ruben Undheim + 2006-2017 Tim Edwards + 2014-2015 Ruben Undheim + 2015 Staf Verhaegen License: GPL-2+_with_Yale_conditions Files: src/Ylib/radixsort.c @@ -18,7 +19,7 @@ Copyright: 1983 Berkeley License: GPL-2+_with_Yale_conditions Files: debian/* -Copyright: 2014-2015 Ruben Undheim +Copyright: 2014-2017 Ruben Undheim License: GPL-2+ diff --git a/debian/libycadgraywolf0.install b/debian/libycadgraywolf1.install similarity index 100% rename from debian/libycadgraywolf0.install rename to debian/libycadgraywolf1.install diff --git a/debian/libycadgraywolf0.symbols b/debian/libycadgraywolf1.symbols similarity index 96% rename from debian/libycadgraywolf0.symbols rename to debian/libycadgraywolf1.symbols index 7e6d42a..cf0d86f 100644 --- a/debian/libycadgraywolf0.symbols +++ b/debian/libycadgraywolf1.symbols @@ -1,4 +1,4 @@ -libycadgraywolf.so.0 libycadgraywolf0 #MINVER# +libycadgraywolf.so.1 libycadgraywolf1 #MINVER# TW3DdrawCube@Base 0.1.1 TW3Dnormal_view@Base 0.1.1 TW3Dperspective@Base 0.1.1 @@ -267,17 +267,6 @@ libycadgraywolf.so.0 libycadgraywolf0 #MINVER# 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 @@ libycadgraywolf.so.0 libycadgraywolf0 #MINVER# 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 --git a/debian/patches/01_no_rpath_in_debian_package.patch b/debian/patches/01_no_rpath_in_debian_package.patch index bb3f58c..4b541f8 100644 --- a/debian/patches/01_no_rpath_in_debian_package.patch +++ b/debian/patches/01_no_rpath_in_debian_package.patch @@ -8,13 +8,13 @@ Forwarded: doesn't make sense upstream 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 --git a/debian/patches/02_dont_create_link.patch b/debian/patches/02_dont_create_link.patch index fdff778..55f1f4e 100644 --- a/debian/patches/02_dont_create_link.patch +++ b/debian/patches/02_dont_create_link.patch @@ -8,11 +8,11 @@ Forwarded: Doesn't make sense upstream 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 --git a/debian/rules b/debian/rules index 253337b..efa65be 100755 --- a/debian/rules +++ b/debian/rules @@ -11,5 +11,5 @@ export DEB_BUILD_MAINT_OPTIONS = hardening=+all dh $@ override_dh_strip: - dh_strip -plibycadgraywolf0 --keep-debug --dbg-package=libycadgraywolf0-dbg + dh_strip -plibycadgraywolf1 --keep-debug --dbg-package=libycadgraywolf1-dbg dh_strip -pgraywolf diff --git a/src/Ylib/CMakeLists.txt b/src/Ylib/CMakeLists.txt index 7e0ce81..f3b9bc0 100644 --- a/src/Ylib/CMakeLists.txt +++ b/src/Ylib/CMakeLists.txt @@ -3,12 +3,12 @@ #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) -set_target_properties(ycadgraywolf PROPERTIES VERSION 0.0.0 SOVERSION 0) +set_target_properties(ycadgraywolf PROPERTIES VERSION 1.0.0 SOVERSION 1) INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/include ${CMAKE_BINARY_DIR}/include) diff --git a/src/Ylib/matrix.c b/src/Ylib/matrix.c deleted file mode 100644 index e9f7a1a..0000000 --- a/src/Ylib/matrix.c +++ /dev/null @@ -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 -#include -#include -#include - -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 --git a/src/Ylib/svd.c b/src/Ylib/svd.c deleted file mode 100644 index c83c50c..0000000 --- a/src/Ylib/svd.c +++ /dev/null @@ -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 -#include -#include -#include - -/* 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 -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 */ diff --git a/src/twmc/CMakeLists.txt b/src/twmc/CMakeLists.txt index 0f4d31b..25bca38 100644 --- a/src/twmc/CMakeLists.txt +++ b/src/twmc/CMakeLists.txt @@ -5,6 +5,7 @@ add_executable(TimberWolfMC acceptt.c analyze.c compact.c config1.c config2.c co 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 --git a/src/twmc/wireratio.c b/src/twmc/wireratio.c index 98b6b62..386a302 100644 --- a/src/twmc/wireratio.c +++ b/src/twmc/wireratio.c @@ -60,10 +60,39 @@ static char SccsId[] = "@(#) wireratio.c version 3.9 3/10/92" ; #include #include #include -#include +#include +#include +#include + +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 @@ INT count ; 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 @@ adapt_wire_estimator() 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 @@ adapt_wire_estimator() 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 @@ adapt_wire_estimator() } 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 @@ adapt_wire_estimator() /* 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 @@ adapt_wire_estimator() 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 @@ adapt_wire_estimator() } /* 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); + } + + /* Compute R = AX - B */ + gsl_vector_memcpy(R, AX); + gsl_vector_sub( R, B ); - /* 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 ) ; + 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 */