Codebase list graywolf / 6d9368c
Merge tag 'upstream/0.1.4+20170306gitecee764' Ruben Undheim 7 years ago
6 changed file(s) with 141 addition(s) and 1040 deletion(s). Raw diff Collapse all Expand all
+8
-0
CMakeLists.txt less more
00 cmake_minimum_required (VERSION 2.6)
11 project (graywolf)
22
3 find_package(PkgConfig)
34 INCLUDE(CheckIncludeFiles)
5
6
7 pkg_check_modules(GSL gsl)
8
9 if (NOT GSL_FOUND)
10 MESSAGE(FATAL_ERROR "The development files for the GNU Scientific Library (libgsl) are required to build graywolf.")
11 endif()
412
513 # Include RPATH in build so that ldconfig is not necessary after install
614 SET(CMAKE_SKIP_BUILD_RPATH FALSE)
+1
-1
src/Ylib/CMakeLists.txt less more
22 #add_executable(mincut main.c output.c readcells.c ${CMAKE_SOURCE_DIR}/src/date/date.c)
33
44
5 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 )
5 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 )
66
77 target_link_libraries(ycadgraywolf X11)
88 target_link_libraries(ycadgraywolf m)
+0
-451
src/Ylib/matrix.c less more
0 /*
1 * Copyright (C) 1991 Yale University
2 *
3 * This work is distributed in the hope that it will be useful; you can
4 * redistribute it and/or modify it under the terms of the
5 * GNU General Public License as published by the Free Software Foundation;
6 * either version 2 of the License,
7 * or any later version, on the following conditions:
8 *
9 * (a) YALE MAKES NO, AND EXPRESSLY DISCLAIMS
10 * ALL, REPRESENTATIONS OR WARRANTIES THAT THE MANUFACTURE, USE, PRACTICE,
11 * SALE OR
12 * OTHER DISPOSAL OF THE SOFTWARE DOES NOT OR WILL NOT INFRINGE UPON ANY
13 * PATENT OR
14 * OTHER RIGHTS NOT VESTED IN YALE.
15 *
16 * (b) YALE MAKES NO, AND EXPRESSLY DISCLAIMS ALL, REPRESENTATIONS AND
17 * WARRANTIES
18 * WHATSOEVER WITH RESPECT TO THE SOFTWARE, EITHER EXPRESS OR IMPLIED,
19 * INCLUDING,
20 * BUT NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
21 * PARTICULAR
22 * PURPOSE.
23 *
24 * (c) LICENSEE SHALL MAKE NO STATEMENTS, REPRESENTATION OR WARRANTIES
25 * WHATSOEVER TO
26 * ANY THIRD PARTIES THAT ARE INCONSISTENT WITH THE DISCLAIMERS BY YALE IN
27 * ARTICLE
28 * (a) AND (b) above.
29 *
30 * (d) IN NO EVENT SHALL YALE, OR ITS TRUSTEES, DIRECTORS, OFFICERS,
31 * EMPLOYEES AND
32 * AFFILIATES BE LIABLE FOR DAMAGES OF ANY KIND, INCLUDING ECONOMIC DAMAGE OR
33 * INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER YALE SHALL BE
34 * ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE
35 * POSSIBILITY OF THE FOREGOING.
36 *
37 */
38
39 /* -----------------------------------------------------------------
40 FILE: matrix.c
41 DESCRIPTION:a curvefitting program.
42 CONTENTS:
43 DATE: Tue Jan 15 01:35:11 EST 1991 - original coding.
44 REVISIONS: Sun Nov 3 12:48:39 EST 1991 - made matrix memory
45 allocation more efficient using YVECTOR routines.
46 12/09/91 - Add prototype macro for non ansi compiler -R.A.Weier
47 ----------------------------------------------------------------- */
48 #ifndef lint
49 static char SccsId[] = "@(#) matrix.c version 1.3 12/9/91" ;
50 #endif
51
52 #include <yalecad/base.h>
53 #include <yalecad/message.h>
54 #include <yalecad/debug.h>
55 #include <yalecad/linalg.h>
56
57 static DOUBLE find_det( P3(DOUBLE **a, INT rows, INT columns ) );
58 static sign_cof( P2(INT row, INT column) ) ;
59
60 /* *****************************************************************
61 Allocate the space for the matrix.
62 **************************************************************** */
63 YMPTR Ymatrix_create( rows, columns )
64 INT rows, columns ;
65 {
66 YMPTR mptr ; /* pointer to matrix record */
67 DOUBLE **m ; /* pointer to matrix space */
68 INT i ;
69
70 mptr = YMALLOC( 1, YMBOX ) ;
71 /* now allocate memory for rows */
72 m = mptr->m = YVECTOR_MALLOC( 1, rows, DOUBLE * ) ;
73 /* now allocate memory for columns */
74 for( i = 1;i<= rows; i++ ){
75 m[i] = YVECTOR_CALLOC( 1, columns, DOUBLE ) ;
76 }
77 mptr->rows = rows ;
78 mptr->columns = columns ;
79 return( mptr ) ;
80 } /* end Ymatrix_create */
81
82 /* *****************************************************************
83 Free the space of a matrix.
84 **************************************************************** */
85 YMPTR Ymatrix_free( mptr )
86 YMPTR mptr ; /* pointer to matrix record */
87 {
88 DOUBLE **m ; /* pointer to matrix space */
89 INT i ;
90 INT rows ;
91
92 rows = mptr->rows ;
93 m = mptr->m ;
94 /* free memory */
95 for( i = 1; i <= rows; i++ ){
96 YVECTOR_FREE( m[i], 1 ) ;
97 }
98 YVECTOR_FREE( m, 1 ) ;
99 YFREE( mptr ) ;
100 } /* end Ymatrix_free */
101
102 YMPTR Ymatrix_transpose( mptr )
103 YMPTR mptr ;
104 {
105 INT i, j ;
106 DOUBLE **m ;
107 DOUBLE **b ;
108 YMPTR buf ;
109
110 buf = Ymatrix_create( mptr->columns, mptr->rows ) ;
111
112 m = mptr->m ;
113 b = buf->m ;
114
115 for( i=1; i <= mptr->rows; i++ ){
116 for( j=1; j <= mptr->columns; j++ ){
117 b[j][i] = m[i][j] ;
118 }
119 }
120 return( buf ) ;
121
122 } /* end Ymatrix_transpose */
123
124 YMPTR Ymatrix_mult( aptr, bptr )
125 YMPTR aptr, bptr ;
126 {
127 INT i, j, k, n ;
128 DOUBLE **a ;
129 DOUBLE **b ;
130 DOUBLE **c ;
131 DOUBLE result ;
132 YMPTR buf ;
133
134 /* first perform error checking */
135 if( aptr->columns != bptr->rows ){
136 sprintf( YmsgG,
137 "Matrices cannot be multiplied together A (%dx%d) B(%dx%d)\n",
138 aptr->rows, aptr->columns, bptr->rows, bptr->columns ) ;
139 M( ERRMSG, "Ymatrix_mult", YmsgG ) ;
140 }
141 buf = Ymatrix_create( aptr->rows, bptr->columns ) ;
142
143 a = aptr->m ;
144 b = bptr->m ;
145 c = buf->m ;
146 n = aptr->columns ;
147
148 for( i=1; i <= aptr->rows; i++ ){
149 for( j=1; j <= bptr->columns; j++ ){
150 result = 0.0 ;
151 for( k=1; k<= n; k++ ){
152 result += a[i][k] * b[k][j] ;
153 }
154 c[i][j] = result ;
155 }
156 }
157 return( buf ) ;
158
159 } /* end Ymatrix_mult */
160
161 YMPTR Ymatrix_sub( aptr, bptr )
162 YMPTR aptr, bptr ;
163 {
164 INT i, j ;
165 DOUBLE **a, **b, **c ;
166 DOUBLE *fast_a, *fast_b, *fast_c ;
167 DOUBLE result ;
168 YMPTR buf ;
169
170 /* first perform error checking */
171 if( aptr->rows != bptr->rows ||
172 aptr->columns != bptr->columns ){
173 sprintf( YmsgG,
174 "Matrices cannot be subtracted A (%dx%d) B(%dx%d)\n",
175 aptr->rows, aptr->columns, bptr->rows, bptr->columns ) ;
176 M( ERRMSG, "Ymatrix_sub", YmsgG ) ;
177 }
178 buf = Ymatrix_create( aptr->rows, aptr->columns ) ;
179
180 a = aptr->m ;
181 b = bptr->m ;
182 c = buf->m ;
183 for( i=1; i <= aptr->rows; i++ ){
184 fast_a = a[i] ;
185 fast_b = b[i] ;
186 fast_c = c[i] ;
187 for( j=1; j <= aptr->columns; j++ ){
188 fast_c[j] = fast_a[j] - fast_b[j] ;
189 }
190 }
191 return( buf ) ;
192
193 } /* end Ymatrix_sub */
194
195 Ymatrix_disp( mptr )
196 YMPTR mptr ;
197 {
198 DOUBLE **m ;
199 INT i, j ;
200 m = mptr->m ;
201
202 for( i=1; i <= mptr->rows; i++ ){
203 for( j=1; j <= mptr->columns; j++ ){
204 fprintf( stderr, "% 4.4le ", m[i][j] ) ;
205 }
206 fprintf( stderr, "\n" ) ;
207 }
208 fprintf( stderr, "\n" ) ;
209 } /* end Ymatrix_disp */
210
211 YMPTR Ymatrix_eye( size )
212 INT size ;
213 {
214 INT i, j ;
215 DOUBLE **m ;
216 DOUBLE *fast ;
217 YMPTR buf ;
218
219 buf = Ymatrix_create( size, size ) ;
220 m = buf->m ;
221 for( i = 1; i <= size; i++ ){
222 fast = m[i] ;
223 for( j = 1; j <= size ; j++ ){
224 if( i == j ){
225 /* one is exact in fp. */
226 fast[j] = 1.0 ;
227 } else {
228 fast[j] = YZERO ;
229 }
230 }
231 }
232 return( buf ) ;
233 } /* end Ymatrix_eye */
234
235 Ymatrix_zero( matrix )
236 YMPTR matrix ;
237 {
238 INT r, c, i, j ;
239 DOUBLE *fast ;
240 DOUBLE **m ;
241
242 m = matrix->m ;
243 r = matrix->rows ;
244 c = matrix->columns ;
245 for( i = 1; i <= r; i++ ){
246 fast = m[i] ;
247 for( j = 1; j <= c ; j++ ){
248 fast[j] = YZERO ;
249 }
250 }
251 } /* end Ymatrix_zero */
252
253
254 YMPTR Ymatrix_copy( input )
255 YMPTR input ;
256 {
257 INT i, j ;
258 DOUBLE **in_mat, **copy_mat ;
259 DOUBLE *in_fast, *copy_fast ;
260 YMPTR copy ;
261
262 copy = Ymatrix_create( input->rows, input->columns ) ;
263 in_mat = input->m ;
264 copy_mat = copy->m ;
265 for( i = 1; i <= input->rows; i++ ){
266 in_fast = in_mat[i] ;
267 copy_fast = copy_mat[i] ;
268 for( j = 1; j <= input->columns ; j++ ){
269 copy_fast[j] = in_fast[j] ;
270 }
271 }
272 return( copy ) ;
273 } /* end Ymatrix_copy */
274
275 YMPTR Ymatrix_linv( aptr )
276 YMPTR aptr ;
277 {
278 INT i, j, k, n ;
279 DOUBLE **a ;
280 DOUBLE **b ;
281 DOUBLE **c ;
282 DOUBLE det, recip_det ;
283 YMPTR cof ;
284 YMPTR buf ;
285
286 buf = Ymatrix_create( aptr->rows, aptr->columns ) ;
287 cof = Ymatrix_cofactors( aptr ) ;
288
289 det = 0.0 ;
290 a = aptr->m ;
291 c = cof->m ;
292 /* find the determinant using first column */
293 for( i=1; i <= aptr->rows; i++ ){
294 det += a[i][1] * c[i][1] ;
295 }
296
297 ASSERT( det == find_det( aptr->m, aptr->rows, aptr->columns ),
298 "Yinv", "Problem with determinant" ) ;
299
300 if( det != 0.0 ){
301 recip_det = 1.0 / det ;
302 buf = Ymatrix_transpose( cof ) ;
303 b = buf->m ;
304 for( i=1; i <= buf->rows; i++ ){
305 for( j=1; j <= buf->columns; j++ ){
306 b[i][j] *= recip_det ;
307 }
308 }
309 }
310 return( buf ) ;
311
312 } /* end Ymatrix_linv */
313
314 YMPTR Ymatrix_cofactors( aptr )
315 YMPTR aptr ;
316 {
317
318 INT i, j, k, l ;
319 INT r, c ;
320 INT rows, columns ;
321 DOUBLE cofactor ;
322 DOUBLE **a, **b, **m ;
323 YMPTR buf ;
324
325 buf = Ymatrix_create( aptr->rows, aptr->columns ) ;
326 rows = aptr->rows ;
327 columns = aptr->columns ;
328 a = aptr->m ; /* a matrix - want cofactors of this matrix */
329 b = buf->m ; /* buf matrix */
330
331 /* allocate space for a row - 1 x column - 1 clone of a matrix */
332 m = YVECTOR_MALLOC( 1, rows - 1, DOUBLE * ) ;
333 /* now allocate memory for columns */
334 for( i = 1;i< rows; i++ ){
335 m[i] = YVECTOR_CALLOC( 1, columns-1, DOUBLE ) ;
336 }
337
338 /* process all cofactors */
339 for( i=1; i <= rows; i++ ){
340 for( j=1; j <= columns; j++ ){
341
342 /* load cofactor array */
343 /* strike row and column */
344 r = 0 ;
345 for( k = 1; k <= rows; k++ ){
346 if( k == i ){
347 continue ;
348 }
349 r++ ;
350 c = 0 ;
351 for( l = 1; l <= columns; l++ ){
352 if( l == j ){
353 continue ;
354 }
355 m[r][++c] = a[k][l] ;
356 }
357 } /* end load of cofactor array */
358
359 /* process one cofactor */
360 cofactor = sign_cof(i,j) * find_det( m, rows-1, columns-1 ) ;
361 /* store result in buf */
362 b[i][j] = cofactor ;
363
364 }
365 }
366 /* now time to free memory */
367 for( i = 1;i< columns; i++ ){
368 YVECTOR_FREE( m[i], 1 ) ;
369 }
370 YVECTOR_FREE( m, 1 ) ;
371 return( buf ) ;
372 } /* end Ymatrix_cofactors */
373
374 /* find a determinant - works on the actual memory arrays - not user */
375 /* records */
376 static DOUBLE find_det( a, rows, columns )
377 DOUBLE **a ;
378 INT rows, columns ;
379 {
380 DOUBLE result ;
381 DOUBLE cofactor ;
382 DOUBLE **m ;
383 INT i, j, k ;
384 INT r, c ;
385
386 if( rows == 1 && columns == 1 ){
387 result = a[1][1] ;
388
389 } else if( rows == 2 && columns == 2 ){
390 result = a[1][1] * a[2][2] - a[2][1] * a[1][2] ;
391
392 } else {
393
394 /* allocate space for an row-1 x columns-1 matrix */
395 m = YVECTOR_MALLOC( 1, rows-1, DOUBLE * ) ;
396 /* now allocate memory for columns */
397 for( i = 1;i< rows; i++ ){
398 m[i] = YVECTOR_MALLOC( 1, columns-1, DOUBLE ) ;
399 }
400
401 /* do Laplace expansion along top row - all columns */
402 result = 0.0 ;
403 for( k=1; k <= columns; k++ ){
404
405 /* load array by striking row and column */
406 /* struck row is always one */
407 r = 0 ;
408 for( i = 2; i <= rows; i++ ){
409 r++ ;
410 c = 0 ;
411 for( j = 1; j <= columns; j++ ){
412 if( j == k ){
413 continue ;
414 }
415 m[r][++c] = a[i][j] ;
416 }
417 } /* end load of cofactor array */
418
419 ASSERT( r == rows-1, "Yfind_det","problem with rows" ) ;
420 ASSERT( c == columns-1, "Yfind_det","problem with rows" ) ;
421
422 /* recursively call find_det to get result */
423 cofactor = sign_cof(1,k) * find_det( m, rows-1, columns-1 ) ;
424 result += a[1][k] * cofactor ;
425 }
426
427 /* now time to free memory */
428 for( i = 1;i< columns; i++ ){
429 YVECTOR_FREE( m[i], 1 ) ;
430 }
431 YVECTOR_FREE( m, 1 ) ;
432 }
433
434 return( result ) ;
435
436 } /* end find_det */
437
438 static sign_cof( row, column )
439 INT row, column ;
440 {
441 INT sum ;
442
443 /* implements (-1) ** (row + column) */
444 sum = row + column ;
445 if( sum % 2 ){ /* odd power */
446 return( -1.0 ) ;
447 } else { /* even power */
448 return( 1.0 ) ;
449 }
450 } /* end sign_cof */
+0
-522
src/Ylib/svd.c less more
0 /*
1 * Copyright (C) 1991 Yale University
2 *
3 * This work is distributed in the hope that it will be useful; you can
4 * redistribute it and/or modify it under the terms of the
5 * GNU General Public License as published by the Free Software Foundation;
6 * either version 2 of the License,
7 * or any later version, on the following conditions:
8 *
9 * (a) YALE MAKES NO, AND EXPRESSLY DISCLAIMS
10 * ALL, REPRESENTATIONS OR WARRANTIES THAT THE MANUFACTURE, USE, PRACTICE,
11 * SALE OR
12 * OTHER DISPOSAL OF THE SOFTWARE DOES NOT OR WILL NOT INFRINGE UPON ANY
13 * PATENT OR
14 * OTHER RIGHTS NOT VESTED IN YALE.
15 *
16 * (b) YALE MAKES NO, AND EXPRESSLY DISCLAIMS ALL, REPRESENTATIONS AND
17 * WARRANTIES
18 * WHATSOEVER WITH RESPECT TO THE SOFTWARE, EITHER EXPRESS OR IMPLIED,
19 * INCLUDING,
20 * BUT NOT LIMITED TO, WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
21 * PARTICULAR
22 * PURPOSE.
23 *
24 * (c) LICENSEE SHALL MAKE NO STATEMENTS, REPRESENTATION OR WARRANTIES
25 * WHATSOEVER TO
26 * ANY THIRD PARTIES THAT ARE INCONSISTENT WITH THE DISCLAIMERS BY YALE IN
27 * ARTICLE
28 * (a) AND (b) above.
29 *
30 * (d) IN NO EVENT SHALL YALE, OR ITS TRUSTEES, DIRECTORS, OFFICERS,
31 * EMPLOYEES AND
32 * AFFILIATES BE LIABLE FOR DAMAGES OF ANY KIND, INCLUDING ECONOMIC DAMAGE OR
33 * INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER YALE SHALL BE
34 * ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE
35 * POSSIBILITY OF THE FOREGOING.
36 *
37 */
38
39 /* -----------------------------------------------------------------
40 FILE: svd.c
41 DESCRIPTION:Singular value decomposition routines.
42 CONTENTS:
43 DATE: Tue Jan 15 00:52:31 EST 1991 - original coding.
44 REVISIONS: Tue Jan 15 22:26:57 PST 1991 - added TRUE return condition.
45 Thu Jan 17 00:59:42 PST 1991 - now notify user of
46 singular value when in debug mode.
47 Thu Jan 24 20:17:23 PST 1991 - changed vector routines.
48 Tue Mar 12 16:56:22 CST 1991 - got rid of unnecessary
49 assignment.
50 ----------------------------------------------------------------- */
51 #ifndef lint
52 static char SccsId[] = "@(#) svd.c version 1.7 12/15/91" ;
53 #endif
54
55 #include <yalecad/base.h>
56 #include <yalecad/message.h>
57 #include <yalecad/debug.h>
58 #include <yalecad/linalg.h>
59
60 /* static definitions */
61 static DOUBLE atS,btS,ctS;
62 static DOUBLE maxarg1S,maxarg2S;
63
64 #undef MAX
65 #define PYTHAG(a,b) ((atS=fabs(a)) > (btS=fabs(b)) ? \
66 (ctS=btS/atS,atS*sqrt(1.0+ctS*ctS)) : (btS ? (ctS=atS/btS,btS*sqrt(1.0+ctS*ctS)): 0.0))
67 #define MAX(a,b) (maxarg1S=(a),maxarg2S=(b),(maxarg1S) > (maxarg2S) ?\
68 (maxarg1S) : (maxarg2S))
69 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
70 #define TOLERANCE 1e-6
71
72 /* ---------------------------------------------------------------
73 Procedure Ysvd_solve - solve a set of linear equations using SVD.
74 Input: A - a m x n matrix.
75 Input: B - a m x k matrix.
76 Output: X - a n x k matrix.
77 -----------------------------------------------------------------*/
78 BOOL Ysvd_solve( A, B, Xret )
79 YMPTR A, B, *Xret ;
80 {
81 INT i, j, jj, k, m, n, t ;
82 YMPTR U, W, V ;
83 DOUBLE wmax, threshold, sum ;
84 DOUBLE **b, **u, **v, **w, **x ;
85 DOUBLE *tmp ;
86
87 m = A->rows ;
88 n = A->columns ;
89 k = B->columns ;
90
91 /* next make sure B conforms to A */
92 if( m != B->rows ){
93 M( ERRMSG, "Ysvd_solve", "The matrices A and B do not conform\n" ) ;
94 sprintf( YmsgG, "A = %d x %d B = %d x %d\n\n", m, n, B->rows, k ) ;
95 M( ERRMSG, NULL, YmsgG ) ;
96 return(FALSE) ;
97 }
98
99 /* first factor A using SVD */
100 if(!(Ysvd_decompose( A, &U, &W, &V ))){
101 return( FALSE ) ;
102 }
103 D( "Ysvd_solve/decompose",
104 YMPTR Vt ;
105 YMPTR WVt ;
106 YMPTR UWVt ;
107
108 fprintf( stderr, "U:\n" ) ;
109 Ymatrix_disp( U ) ;
110 fprintf( stderr, "W:\n" ) ;
111 Ymatrix_disp( W ) ;
112 fprintf( stderr, "V:\n" ) ;
113 Ymatrix_disp( V ) ;
114
115 Vt = Ymatrix_transpose( V ) ;
116 fprintf( stderr, "Vt:\n" ) ;
117 Ymatrix_disp( Vt ) ;
118 WVt = Ymatrix_mult( W, Vt ) ;
119 UWVt = Ymatrix_mult( U, WVt ) ;
120 fprintf( stderr, "U * W * Vt:\n" ) ;
121 Ymatrix_disp( UWVt ) ;
122 ) ;
123
124 /* now throw away those things which make the answer unstable */
125 /* first find the maximum w */
126 wmax=0.0;
127 w = W->m ;
128 for (j=1;j<=n;j++){
129 if(w[j][j] > wmax){
130 wmax=w[j][j] ;
131 }
132 }
133 /* next create a threshold */
134 threshold = TOLERANCE * wmax ;
135 /* now zero below threshold */
136 for (j=1;j<=n;j++){
137 if( w[j][j] < threshold ){
138 D( "Ysvd_solve/singular",
139 fprintf( stderr, "Singular value for w[%d]:%4.6le\n",
140 j, w[j][j] ) ;
141 ) ;
142 w[j][j] = 0.0;
143 }
144 }
145 D( "Ysvd_solve/a_weight",
146 fprintf( stderr, "W:\n" ) ;
147 Ymatrix_disp( W ) ;
148 ) ;
149
150 D( "Ysvd_solve/backsub",
151 YMPTR Ut ;
152 YMPTR Utb ;
153 YMPTR WUtb ;
154 YMPTR Xver ;
155 YMPTR Wprime ;
156 DOUBLE **wp ;
157
158 Ut = Ymatrix_transpose( U ) ;
159 Utb = Ymatrix_mult( Ut, B ) ;
160 Wprime = Ymatrix_copy( W ) ;
161 wp = Wprime->m ;
162 for( i = 1; i <= n; i++ ){
163 if( wp[i][i] ){
164 wp[i][i] = 1 / wp[i][i] ;
165 }
166 }
167 WUtb = Ymatrix_mult( Wprime, Utb ) ;
168 Xver = Ymatrix_mult( V, WUtb ) ;
169 fprintf( stderr, "V * 1 / w * Ut * B:\n" ) ;
170 Ymatrix_disp( Xver ) ;
171 ) ;
172
173 /* now we are ready for back substitution */
174 /* x = V * diag( 1 / w ) * Ut * b */
175 *Xret = Ymatrix_create( n, k ) ;
176 x = (*Xret)->m ;
177 u = U->m ;
178 v = V->m ;
179 b = B->m ;
180 tmp = YVECTOR_MALLOC(1,n,DOUBLE) ;
181 for( t = 1; t <= k ; t++ ){
182 /* first calculate Ut * b */
183 for (j=1;j<=n;j++) {
184 sum=0.0;
185 if (w[j][j]) {
186 for (i=1;i<=m;i++){
187 sum += u[i][j]*b[i][t];
188 }
189 /* this is the divide by 1 / w */
190 sum /= w[j][j];
191 }
192 tmp[j]=sum;
193 }
194 /* now multiply by V */
195 for (j=1;j<=n;j++) {
196 sum=0.0;
197 for (jj=1;jj<=n;jj++){
198 sum += v[j][jj]*tmp[jj];
199 }
200 x[j][t]=sum;
201 }
202 } /* end looping on columns */
203
204 /* free memory */
205 YVECTOR_FREE(tmp,1) ;
206 Ymatrix_free( U ) ;
207 Ymatrix_free( V ) ;
208 Ymatrix_free( W ) ;
209 return( TRUE ) ;
210
211 } /* end Ysvd_solve */
212
213 /* -------------------------------------------------------------
214 Procedure Ysvd_deccompose - decompose A into U W Vt.
215 Input: A an m x n matrix.
216 Output: U - a m x n column orthogonal matrix.
217 Output: W - a n x 1 column vector.
218 Output: V - a n x n column orthogonal matrix.
219 ---------------------------------------------------------------- */
220 BOOL Ysvd_decompose( A, Uret, Wret, Vret )
221 YMPTR A, *Uret, *Wret, *Vret ;
222 {
223 INT m, n ;
224 INT flag,i,its,j,jj,k,l,nm;
225 DOUBLE c,f,h,s,x,y,z;
226 DOUBLE anorm=0.0,g=0.0,scale=0.0;
227 DOUBLE *rv1;
228 DOUBLE **a,*w,**v, **wfast ;
229
230 m = A->rows ;
231 n = A->columns ;
232
233 if (m < n){
234 sprintf( YmsgG, "m < n (A = %d x %d)\n", m, n ) ;
235 M( ERRMSG, "Ysvd_decompose", YmsgG ) ;
236 M( ERRMSG, NULL, "You must augment A with extra zero rows");
237 return( FALSE ) ;
238 }
239
240 /* create the working space */
241 *Uret = Ymatrix_copy( A ) ;
242 a = (*Uret)->m ;
243 *Wret = Ymatrix_eye( n ) ;
244 *Vret = Ymatrix_create( n, n ) ;
245 v = (*Vret)->m ;
246 w = YVECTOR_MALLOC(1,n,DOUBLE);
247 rv1 = YVECTOR_MALLOC(1,n,DOUBLE);
248
249 /* Householder reduction to bidiagonal form */
250 for (i=1;i<=n;i++) {
251 l=i+1;
252 rv1[i]=scale*g;
253 g=s=scale=0.0;
254 if (i <= m) {
255 for (k=i;k<=m;k++) scale += fabs(a[k][i]);
256 if (scale) {
257 for (k=i;k<=m;k++) {
258 a[k][i] /= scale;
259 s += a[k][i]*a[k][i];
260 }
261 f=a[i][i];
262 g = -SIGN(sqrt(s),f);
263 h=f*g-s;
264 a[i][i]=f-g;
265 if (i != n) {
266 for (j=l;j<=n;j++) {
267 for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
268 f=s/h;
269 for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
270 }
271 }
272 for (k=i;k<=m;k++) a[k][i] *= scale;
273 }
274 }
275 w[i]=scale*g;
276 g=s=scale=0.0;
277 if (i <= m && i != n) {
278 for (k=l;k<=n;k++) scale += fabs(a[i][k]);
279 if (scale) {
280 for (k=l;k<=n;k++) {
281 a[i][k] /= scale;
282 s += a[i][k]*a[i][k];
283 }
284 f=a[i][l];
285 g = -SIGN(sqrt(s),f);
286 h=f*g-s;
287 a[i][l]=f-g;
288 for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
289 if (i != m) {
290 for (j=l;j<=m;j++) {
291 for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
292 for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
293 }
294 }
295 for (k=l;k<=n;k++) a[i][k] *= scale;
296 }
297 }
298 anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i])));
299 } /* end of Householder reduction */
300
301 /* Accummulation of right-hand transformations */
302 for (i=n;i>=1;i--) {
303 if (i < n) {
304 if (g) {
305 for (j=l;j<=n;j++)
306 v[j][i]=(a[i][j]/a[i][l])/g;
307 for (j=l;j<=n;j++) {
308 for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
309 for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
310 }
311 }
312 for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
313 }
314 v[i][i]=1.0;
315 g=rv1[i];
316 l=i;
317 } /* end accummulation of right-hand transformation */
318
319 /* Accummulation of left-hand transformations */
320 for (i=n;i>=1;i--) {
321 l=i+1;
322 g=w[i];
323 if (i < n) for (j=l;j<=n;j++) a[i][j]=0.0;
324 if (g) {
325 g=1.0/g;
326 if (i != n) {
327 for (j=l;j<=n;j++) {
328 for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
329 f=(s/a[i][i])*g;
330 for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
331 }
332 }
333 for (j=i;j<=m;j++) a[j][i] *= g;
334 } else {
335 for (j=i;j<=m;j++) a[j][i]=0.0;
336 }
337 ++a[i][i];
338 } /* end accummulation of left-hand transformation */
339
340 /* Diagonalization of the bidiagonal form */
341 /* loop over singular values */
342 for (k=n;k>=1;k--) {
343 /* loop over allowed iterations */
344 for (its=1;its<=30;its++) {
345 flag=1;
346 for (l=k;l>=1;l--) {
347 nm=l-1;
348 if (fabs(rv1[l])+anorm == anorm) {
349 flag=0;
350 break;
351 }
352 if (fabs(w[nm])+anorm == anorm) break;
353 }
354 if (flag) {
355 s=1.0;
356 for (i=l;i<=k;i++) {
357 f=s*rv1[i];
358 if (fabs(f)+anorm != anorm) {
359 g=w[i];
360 h=PYTHAG(f,g);
361 w[i]=h;
362 h=1.0/h;
363 c=g*h;
364 s=(-f*h);
365 for (j=1;j<=m;j++) {
366 y=a[j][nm];
367 z=a[j][i];
368 a[j][nm]=y*c+z*s;
369 a[j][i]=z*c-y*s;
370 }
371 }
372 }
373 }
374 /* Convergence */
375 z=w[k];
376 if (l == k) {
377 if (z < 0.0) {
378 w[k] = -z;
379 for (j=1;j<=n;j++) v[j][k]=(-v[j][k]);
380 }
381 break;
382 }
383 if (its == 30){
384 M(ERRMSG, "Ysvd_decompose", "No convergence in 30 SVDCMP iterations");
385 return( FALSE ) ;
386 }
387 x=w[l];
388 nm=k-1;
389 y=w[nm];
390 g=rv1[nm];
391 h=rv1[k];
392 f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
393 g=PYTHAG(f,1.0);
394 f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
395
396 /* The QR transformation */
397 c=s=1.0;
398 for (j=l;j<=nm;j++) {
399 i=j+1;
400 g=rv1[i];
401 y=w[i];
402 h=s*g;
403 g=c*g;
404 z=PYTHAG(f,h);
405 rv1[j]=z;
406 c=f/z;
407 s=h/z;
408 f=x*c+g*s;
409 g=g*c-x*s;
410 h=y*s;
411 y=y*c;
412 for (jj=1;jj<=n;jj++) {
413 x=v[jj][j];
414 z=v[jj][i];
415 v[jj][j]=x*c+z*s;
416 v[jj][i]=z*c-x*s;
417 }
418 z=PYTHAG(f,h);
419 w[j]=z;
420 if (z) {
421 z=1.0/z;
422 c=f*z;
423 s=h*z;
424 }
425 f=(c*g)+(s*y);
426 x=(c*y)-(s*g);
427 for (jj=1;jj<=m;jj++) {
428 y=a[jj][j];
429 z=a[jj][i];
430 a[jj][j]=y*c+z*s;
431 a[jj][i]=z*c-y*s;
432 }
433 }
434 rv1[l]=0.0;
435 rv1[k]=f;
436 w[k]=x;
437 }
438 }
439
440 /* now transfer w to Wret */
441 wfast = (*Wret)->m ;
442 for( i = 1; i <= n; i++ ){
443 wfast[i][i] = w[i] ;
444 }
445 YVECTOR_FREE(rv1,1) ;
446 YVECTOR_FREE(w,1) ;
447 return( TRUE ) ;
448
449 } /* end Ysvd_decompose */
450
451 #undef SIGN
452 #undef MAX
453 #undef PYTHAG
454
455 #ifdef TEST
456 /* *************************** TEST SOLVER ****************** */
457 #include <yalecad/cleanup.h>
458 main(argc,argv)
459 int argc;
460 char *argv[];
461 {
462
463 YMPTR X, y, Xt, XtX, B, Xret ;
464
465 /* start up cleanup handler */
466 YINITCLEANUP( argv[0], NULL, MAYBEDUMP ) ;
467 YsetDebug( TRUE ) ;
468
469 X = Ymatrix_create( 5, 3 ) ;
470 y = Ymatrix_create( 5, 1 ) ;
471
472 /* load matrix */
473 YMATRIX( X, 1, 1 ) = 1.0 ;
474 YMATRIX( X, 1, 2 ) = 0.0 ;
475 YMATRIX( X, 1, 3 ) = 0.0 ;
476 YMATRIX( X, 2, 1 ) = 1.0 ;
477 YMATRIX( X, 2, 2 ) = 67.0 ;
478 YMATRIX( X, 2, 3 ) = 448.0 ;
479 YMATRIX( X, 3, 1 ) = 1.0 ;
480 YMATRIX( X, 3, 2 ) = 157.0 ;
481 YMATRIX( X, 3, 3 ) = 24649.0 ;
482 YMATRIX( X, 4, 1 ) = 1.0 ;
483 YMATRIX( X, 4, 2 ) = 67.0 ;
484 YMATRIX( X, 4, 3 ) = 448.0 ;
485 YMATRIX( X, 5, 1 ) = 1.0 ;
486 YMATRIX( X, 5, 2 ) = 0.0 ;
487 YMATRIX( X, 5, 3 ) = 0.0 ;
488 M( MSG, NULL, "X:\n" ) ;
489 Ymatrix_disp( X ) ;
490
491 YMATRIX( y, 1, 1 ) = 0.7 ;
492 YMATRIX( y, 2, 1 ) = 0.7 ;
493 YMATRIX( y, 3, 1 ) = 0.15 ;
494 YMATRIX( y, 4, 1 ) = 0.7 ;
495 YMATRIX( y, 5, 1 ) = 0.7 ;
496 M( MSG, NULL, "y:\n" ) ;
497 Ymatrix_disp( y ) ;
498
499 M( MSG, NULL, "\nXt:\n" ) ;
500 Xt = Ymatrix_transpose( X ) ;
501 Ymatrix_disp( Xt ) ;
502
503 M( MSG, NULL, "\nXtX:\n" ) ;
504 XtX = Ymatrix_mult( Xt, X ) ;
505 Ymatrix_disp( XtX ) ;
506
507 /* now test eq solver */
508 B = Ymatrix_mult( Xt, y ) ;
509 M( MSG, NULL, "\nB originally\n");
510 Ymatrix_disp( B ) ;
511
512 if( Ysvd_solve( XtX, B, &Xret )){
513 M( MSG, NULL, "\nSvd Solver\n");
514 Ymatrix_disp( Xret ) ;
515 }
516
517 YexitPgm(0) ;
518
519 } /* end main program */
520
521 #endif /* TEST */
+1
-0
src/twmc/CMakeLists.txt less more
44 target_link_libraries(TimberWolfMC ${CMAKE_BINARY_DIR}/src/Ylib/libycadgraywolf.so)
55 target_link_libraries(TimberWolfMC X11)
66 target_link_libraries(TimberWolfMC m)
7 target_link_libraries(TimberWolfMC ${GSL_LIBRARIES})
78
89 INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/include ${CMAKE_BINARY_DIR}/include .)
910
+131
-66
src/twmc/wireratio.c less more
5959 #include <dens.h>
6060 #include <yalecad/debug.h>
6161 #include <yalecad/file.h>
62 #include <yalecad/linalg.h>
62 #include <gsl/gsl_vector.h>
63 #include <gsl/gsl_matrix.h>
64 #include <gsl/gsl_linalg.h>
65
66 gsl_matrix_disp( mptr, rows, cols )
67 gsl_matrix *mptr ;
68 int rows, cols;
69 {
70 INT i, j ;
71
72 for( i=0; i < rows; i++ ){
73 for( j=0; j < cols; j++ ){
74 fprintf( stderr, "% 4.4le ", gsl_matrix_get(mptr, i, j)) ;
75 }
76 fprintf( stderr, "\n" ) ;
77 }
78 fprintf( stderr, "\n" ) ;
79 } /* end gsl_matrix_disp */
80
81 gsl_vector_disp( vptr, rows )
82 gsl_vector *vptr ;
83 int rows;
84 {
85 INT i;
86
87 for( i=0; i < rows; i++ ){
88 fprintf( stderr, "% 4.4le ", gsl_vector_get(vptr, i)) ;
89 }
90 fprintf( stderr, "\n" ) ;
91 } /* end gsl_vector_disp */
6392
6493 static set_pins( A, center, loc, tile_side, sidepins, count )
65 YMPTR A ; /* the matrix holding x y positions */
94 gsl_matrix *A ; /* the matrix holding x y positions */
6695 INT center ;
6796 INT loc ;
6897 INT tile_side ;
74103 if( sidepins ){
75104 side = find_tile_side( center, loc, tile_side ) ;
76105 if( side ){
77 YMATRIX( A, count, 6 ) = (DOUBLE) sidepins[side] ;
106 gsl_matrix_set(A, count, 5, (DOUBLE) sidepins[side]);
78107 } else {
79108 M( ERRMSG, "adapt_wire_estimator",
80109 "Trouble finding pinside - defaulting to 0.\n" ) ;
81 YMATRIX( A, count, 6 ) = 0.0 ;
110 gsl_matrix_set(A, count, 5, 0.0);
82111 }
83112 } else {
84113 /* no pins for cell */
85 YMATRIX( A, count, 6 ) = 0.0 ;
114 gsl_matrix_set(A, count, 5, 0.0);
86115 }
87116 } /* end set_pins */
88117
94123 INT xc, yc ; /* cell center */
95124 INT *sidepins ; /* array holding #pins for side */
96125 INT l, r, b, t ; /* the global position of the rtiles */
126 INT solved ; /* status of gsl_linalg_SV_solve */
97127 INT *find_pin_sides() ; /* find number of pins on all sides */
98128 char filename[LRECL] ; /* output the results of the SVD fit */
99129 FILE *fp ; /* write out the results */
100 YMPTR A ; /* the matrix holding x y positions */
101 YMPTR B ; /* the resulting global routing space*/
102 YMPTR Xret ; /* the solution to At*A*x = At*B */
130 gsl_matrix *A ; /* the matrix holding x y positions */
131 gsl_vector *B ; /* the resulting global routing space*/
132 gsl_vector *Xret ; /* the solution to At*A*x = At*B */
103133 DOUBLE x, y ; /* cell placement */
104134 DOUBLE lf, rf, bf, tf ; /* the global position of the rtiles */
105135 DOUBLE xlength, ylength;/* length of side */
106136 DOUBLE xrouting, yrouting;/* routing space */
107137 CELLBOXPTR cptr ; /* current pointer to cell */
108138 RTILEBOXPTR rptr ; /* traverse tiles */
139
140 gsl_matrix *U ;
141 gsl_matrix *V ;
142 gsl_vector *S ;
143 gsl_vector *work ;
109144
110145 if(!(routingTilesG)){
111146 return ;
123158 }
124159 if( count < 6 ){
125160 /* now we can size the matrices */
126 A = Ymatrix_create( 6, 6 ) ;
127 B = Ymatrix_create( 6, 1 ) ;
161 A = gsl_matrix_alloc(6, 6);
162 B = gsl_vector_alloc(6);
128163 } else {
129164 /* now we can size the matrices */
130 A = Ymatrix_create( count, 6 ) ;
131 B = Ymatrix_create( count, 1 ) ;
165 A = gsl_matrix_alloc(count, 6);
166 B = gsl_vector_alloc(count);
132167 }
133168
134169 /* initialize the pin counting routines */
150185 /* switchbox densities are not accurate discard */
151186 continue ;
152187 }
153 count++ ;
154 YMATRIX( A, count, 1 ) = 1.0 ; /* for finding const */
188 gsl_matrix_set(A, count, 0, 1.0); /* for finding const */
155189 l = xc + rptr->x1 ;
156190 r = xc + rptr->x2 ;
157191 b = yc + rptr->y1 ;
170204 x = rf ;
171205 y = (bf + tf) / 2.0 ;
172206 set_pins( A, (b+t)/2, r, TILEL, sidepins, count ) ;
173 YMATRIX( A, count, 2 ) = x ;
174 YMATRIX( A, count, 3 ) = x * x ;
175 YMATRIX( A, count, 4 ) = y ;
176 YMATRIX( A, count, 5 ) = y * y ;
177 YMATRIX( B, count, 1 ) = xrouting ;
207 gsl_matrix_set(A, count, 1, x);
208 gsl_matrix_set(A, count, 2, x * x);
209 gsl_matrix_set(A, count, 3, y);
210 gsl_matrix_set(A, count, 4, y * y);
211 gsl_vector_set(B, count, xrouting);
178212 break ;
179213 case TILET:
180214 /* calculate x and y */
181215 x = (lf + rf) / 2.0 ;
182216 y = bf ;
183217 set_pins( A, (l+r)/2, b, TILET, sidepins, count ) ;
184 YMATRIX( A, count, 2 ) = x ;
185 YMATRIX( A, count, 3 ) = x * x ;
186 YMATRIX( A, count, 4 ) = y ;
187 YMATRIX( A, count, 5 ) = y * y ;
188 YMATRIX( B, count, 1 ) = yrouting ;
218 gsl_matrix_set(A, count, 1, x);
219 gsl_matrix_set(A, count, 2, x * x);
220 gsl_matrix_set(A, count, 3, y);
221 gsl_matrix_set(A, count, 4, y * y);
222 gsl_vector_set(B, count, yrouting);
189223 break ;
190224 case TILER:
191225 /* calculate x and y */
192226 x = lf ;
193227 y = (bf + tf) / 2.0 ;
194228 set_pins( A, (b+t)/2, l, TILER, sidepins, count ) ;
195 YMATRIX( A, count, 2 ) = x ;
196 YMATRIX( A, count, 3 ) = x * x ;
197 YMATRIX( A, count, 4 ) = y ;
198 YMATRIX( A, count, 5 ) = y * y ;
199 YMATRIX( B, count, 1 ) = xrouting ;
229 gsl_matrix_set(A, count, 1, x);
230 gsl_matrix_set(A, count, 2, x * x);
231 gsl_matrix_set(A, count, 3, y);
232 gsl_matrix_set(A, count, 4, y * y);
233 gsl_vector_set(B, count, xrouting);
200234 break ;
201235 case TILEB:
202236 /* calculate x and y */
203237 x = (lf + rf) / 2.0 ;
204238 y = tf ;
205239 set_pins( A, (l+r)/2, t, TILEB, sidepins, count ) ;
206 YMATRIX( A, count, 2 ) = x ;
207 YMATRIX( A, count, 3 ) = x * x ;
208 YMATRIX( A, count, 4 ) = y ;
209 YMATRIX( A, count, 5 ) = y * y ;
210 YMATRIX( B, count, 1 ) = yrouting ;
240 gsl_matrix_set(A, count, 1, x);
241 gsl_matrix_set(A, count, 2, x * x);
242 gsl_matrix_set(A, count, 3, y);
243 gsl_matrix_set(A, count, 4, y * y);
244 gsl_vector_set(B, count, yrouting);
211245 break ;
212246 } /* end switch */
247 count++ ;
213248 }
214249 if( sidepins ){
215250 Yvector_free( sidepins, 1, sizeof(INT) ) ;
218253 } /* end loop on cells */
219254
220255 /* now make sure we have at least 6 rows */
221 if( count <= 1 ){
256 if( count < 1 ){
222257 M( ERRMSG, "adapt_wire_estimator",
223258 "Too few cells for TimberWolfMC. Must abort\n" ) ;
224259 YexitPgm(PGMFAIL) ;
225260
226261 } else if( count < 6 ){
227262 /* pad with zeros */
228 for( ++count; count <= 6; count++ ){
229 for( i = 1; i <= 6; i++ ){
230 YMATRIX( A, count, i ) = YMATRIX( A, 1, i ) ;
263 for( ; count < 6; count++ ){
264 for( i = 0; i < 6; i++ ){
265 gsl_matrix_set(A, count, i, gsl_matrix_get(A, 0, i));
231266 }
232 YMATRIX( B, count, 1 ) = YMATRIX( B, 1, 1 ) ;
267 gsl_vector_set(B, count, gsl_vector_get(B, 0));
233268 }
234269 }
235270
236271 D( "TWMC/awe/init",
237272 fprintf( stderr,"\nA:\n" ) ;
238 Ymatrix_disp( A ) ;
273 gsl_matrix_disp( A, count < 6 ? 6 : count, 6 ) ;
239274 fprintf( stderr,"\nB:\n" ) ;
240 Ymatrix_disp( B ) ;
275 gsl_vector_disp( B, count < 6 ? 6 : count ) ;
241276 ) ;
242277
243278 /* now solve using SVD */
244 if( Ysvd_solve( A, B, &Xret )){
279 U = gsl_matrix_alloc((count < 6) ? 6 : count, 6);
280 V = gsl_matrix_alloc(6, 6);
281 S = gsl_vector_alloc(6);
282 work = gsl_vector_alloc(6);
283 gsl_matrix_memcpy(U, A);
284 gsl_linalg_SV_decomp(U, V, S, work);
285 solved = gsl_linalg_SV_solve(U, V, S, B, Xret);
286
287 if( solved ){
245288 sprintf( filename, "%s.mest", cktNameG ) ;
246289 fp = TWOPEN( filename, "w", ABORT ) ;
247290 fprintf( fpoG, "\nThe results of the SVD fit are:\n" ) ;
248 for( i = 1; i <= 6; i++ ){
249 HPO( fp, YMATRIX(Xret,i,1) ) ;
250 HPO( fpoG, YMATRIX(Xret,i,1) ) ;
291 for( i = 0; i < 6; i++ ){
292 HPO( fp, gsl_vector_get(Xret, i));
293 HPO( fpoG, gsl_vector_get(Xret, i));
251294 }
252295 TWCLOSE( fp ) ;
253296 }
254297
298 gsl_matrix_free(U);
299 gsl_matrix_free(V);
300 gsl_vector_free(S);
301 gsl_vector_free(work);
255302
256303 D( "TWMC/awe/answer",
257 YMPTR AX ; /* multiply the answer */
258 YMPTR R ;
259 YMPTR Q ;
304 gsl_vector *AX ; /* multiply the answer */
305 gsl_vector *R ;
306 gsl_matrix *Q ;
260307 INT c ;
261
262 /* now output the differences */
263 AX = Ymatrix_mult( A, Xret ) ;
264 R = Ymatrix_sub( AX, B ) ;
265 Q = Ymatrix_create( count, 7 ) ;
266 for( c = 1; c <= count; c++ ){
267 YMATRIX( Q, c, 1 ) = YMATRIX( A, c, 2 ) ; /* x */
268 YMATRIX( Q, c, 2 ) = YMATRIX( A, c, 4 ) ; /* y */
269 YMATRIX( Q, c, 3 ) = YMATRIX( A, c, 6 ) ; /* pins */
270 YMATRIX( Q, c, 4 ) = YMATRIX( B, c, 1 ) ; /* B */
271 YMATRIX( Q, c, 5 ) = YMATRIX( AX, c, 1 ) ; /* AX */
272 YMATRIX( Q, c, 6 ) = YMATRIX( R, c, 1 ) ; /* B - AX */
273 YMATRIX( Q, c, 7 ) = /* error */
274 YMATRIX( Q, c, 6 ) / YMATRIX( Q, c, 4 ) ;
308 DOUBLE d ;
309 DOUBLE dx ;
310
311 AX = gsl_vector_alloc(count < 6 ? count : 6);
312 R = gsl_vector_alloc(count < 6 ? count : 6);
313
314 /* Compute AX = A * Xret */
315 for ( c = 0; c < (count < 6) ? 6 : count; c++ ) {
316 d = 0.0;
317 dx = gsl_vector_get(Xret, c);
318 for ( i = 0; i < 6; i++ ) {
319 d += gsl_matrix_get(A, c, i) * dx;
320 gsl_set_vector(AX, c, d);
321 }
322
323 /* Compute R = AX - B */
324 gsl_vector_memcpy(R, AX);
325 gsl_vector_sub( R, B );
326
327 Q = gsl_matrix_alloc( count, 7 );
328 for( c = 0; c < count; c++ ){
329 gsl_matrix_set(Q, c, 0, gsl_matrix_get(A, c, 1)) ; /* x */
330 gsl_matrix_set(Q, c, 1, gsl_matrix_get(A, c, 3)) ; /* y */
331 gsl_matrix_set(Q, c, 2, gsl_matrix_get(A, c, 5)) ; /* pins */
332 gsl_matrix_set(Q, c, 3, gsl_vector_get(B, c)) ; /* B */
333 gsl_matrix_set(Q, c, 4, gsl_matrix_get(AX, c, 0)) ; /* AX */
334 gsl_matrix_set(Q, c, 5, gsl_matrix_get(R, c, 0)) ; /* B - AX */
335 gsl_matrix_set(Q, c, 6, gsl_matrix_get(Q, c, 5) /
336 gsl_matrix_get(Q, c, 3)) ; /* error */
275337 }
276338 fprintf( stderr, "\n x, y, B, pins, AB, B - AX, error:\n");
277 Ymatrix_disp( Q ) ;
339 gsl_matrix_disp( Q, count, 7 ) ;
340 gsl_matrix_free(Q);
341 gsl_vector_free(R);
342 gsl_vector_free(AX);
278343 ) ;
279344
280345 /* now done free the matrices */
281 Ymatrix_free( A ) ;
282 Ymatrix_free( B ) ;
283 Ymatrix_free( Xret ) ;
346 gsl_matrix_free( A );
347 gsl_vector_free( B );
348 gsl_vector_free( Xret );
284349
285350 } /* end adapt_wire_estimator */