mm-260 === Subject: Where can I find GSL examples? Hi! I'm engineer and amateur mathematician. I'm usually used Gnu Scientific Libray. But GSL tutorial examples is so-so! Where can I find GSL example for Engineering Mathematics and physics? I should be most grateful if you could this, or for any suggestions you may care to take. === Subject: Re: Optimisation where constraint A OR constraint B must be true > I have a convex optimisation problem with a quadratic objective, a few > linear constraints and a lot of quadratic constraints. I can solve this > using SQP (specifically NAG routine E04UCF, which I think is the same as > NPSOL). > > I want to extend this such that each of the quadratic contraints is replaced > by two quadratic constraints. However, I do not need both constraints to be > satisfied, it is suffiecient if either one or the other (or both) is/are > satisfied. > This is called a disjunctive constraint; maybe this helps you search > for it. Disjunctive constraints destroy the convexity of your > problem; the standard way to model them is to introduce additional > binary variables. > Can anyone point me to suitable software or literature? > Maybe my survey > Complete Search in Continuous Global Optimization and > Constraint Satisfaction, > http://www.mat.univie.ac.at/~neum/papers.html#glopt03 > helps. That is very helpful, thank you. As you said, just knowing it is called a disjunctive constraint helps. Just to clarify, by introducing additional binary variables do you mean something like the following: If the constraint is: A(x)>=0 OR B(x)>=0 I introduce two binary variables bA and bB and use constraints: bA*A(x)>=0 AND bB*B(x)>=0 AND bA+bB>=1 -- Dr Andrew McLean Marine and Acoustics Centre QinetiQ ltd The views expressed above are entirely those of the writer and do not represent the views, policy or understanding of any other person or official body. === Subject: Re: Optimisation where constraint A OR constraint B must be true bA*A(x)>=0 bB*B(x)>=0 bA+bB>=1 bA,bB in {0,1} is not always the best formulation. One problematic issue is that when bA=0, A(x) vanishes. Such constructs are sometimes causing difficulties in NLP solvers. An alternative formulation could be: A(x)>= M1*bA B(x)>= M2*(1-bA) bA in {0,1} where M1 = lowerbound on A(x) and M2 = lowerbound on B(x) This formulation does not increase the nonlinearity of the problem. Of course it remains an MINLP. -------------------------------------------------------------- -- Erwin Kalvelagen erwin@gams.com, http://www.gams.com/~erwin -------------------------------------------------------------- -- >> I have a convex optimisation problem with a quadratic objective, a few >> linear constraints and a lot of quadratic constraints. I can solve this >> using SQP (specifically NAG routine E04UCF, which I think is the same as >> NPSOL). >> >> I want to extend this such that each of the quadratic contraints is > replaced >> by two quadratic constraints. However, I do not need both constraints to > be >> satisfied, it is suffiecient if either one or the other (or both) is/are >> satisfied. >> This is called a disjunctive constraint; maybe this helps you search >> for it. Disjunctive constraints destroy the convexity of your >> problem; the standard way to model them is to introduce additional >> binary variables. >> Can anyone point me to suitable software or literature? >> Maybe my survey >> Complete Search in Continuous Global Optimization and >> Constraint Satisfaction, >> http://www.mat.univie.ac.at/~neum/papers.html#glopt03 >> helps. > That is very helpful, thank you. As you said, just knowing it is called a > disjunctive constraint helps. > Just to clarify, by introducing additional binary variables do you mean > something like the following: > If the constraint is: > A(x)>=0 OR B(x)>=0 > I introduce two binary variables bA and bB and use constraints: > bA*A(x)>=0 AND bB*B(x)>=0 AND bA+bB>=1 > -- > Dr Andrew McLean > Marine and Acoustics Centre > QinetiQ ltd > The views expressed above are entirely those of the writer and do > not represent the views, policy or understanding of any other > person or official body. === Subject: Re: Optimisation where constraint A OR constraint B must be true >That is very helpful, thank you. As you said, just knowing it is called a >disjunctive constraint helps. > >Just to clarify, by introducing additional binary variables do you mean >something like the following: > >If the constraint is: > >A(x)>=0 OR B(x)>=0 > >I introduce two binary variables bA and bB and use constraints: > >bA*A(x)>=0 AND bB*B(x)>=0 AND bA+bB>=1 > and bA in {0,1} and bB in {0,1} binary variables this makes the problem to a MINLP mixed integer nonlinear program which is usually solved by branch and bound , in the manner global optimization problems are solved (-> BARON, COCONUT, ..... ) you might get the idea to use ba*(1-bA)=0 and bB*(1-bB)=0 but this does not solve it since E04UCF and others like this one give a local solution only, this is as if you had posed just randomly one of your constraints hth peter === Subject: Re: Optimisation where constraint A OR constraint B must be true > Just to clarify, by introducing additional binary variables do you mean > something like the following: > > If the constraint is: > > A(x)>=0 OR B(x)>=0 > > I introduce two binary variables bA and bB and use constraints: > > bA*A(x)>=0 AND bB*B(x)>=0 AND bA+bB>=1 > Yes; that's one way to do it. There is an extended literature on disjunctive programming, and the naive reformulation is not always the best one. So I'd recommend you do some reading. Arnold Neumaier === Subject: Vortex isolation subroutine by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i29CwFS06382; I have a large data and I want to isolate one vortex from many vortices inside a box. If there any subroutine help me to isolate one vortex, it will be helpful to me. I hope that I can isolate and track vorticity features at different time steps. === Subject: Re: Examples using IMSL for non-linear curve fitting > > The IMSL for FORTRAN library has routines to perform non-linear least square > curve fitting. Is there any example FORTRAN code using IMSL for non-linear > curve fitting available so I can learn how to use the library. I will use > CVF 6.6. > Unsolicited advice: beware of using the IMSL, NAG or other proprietary libraries for algorithms that have also been implemented well in the public domain. Non-linear least squares falls into this category -- codes for this have been discussed several times on sci.math.num-analysis. Compaq Visual Fortran is a very good compiler, but if you use IMSL libraries of the professional vesrsion, you are locked into CVF and cannot test your code with other compileUNLESS you are willing to pay several hundred dollars for the IMSL license for each new compiler. And there are compilers for which IMSL is just not available. The successor to CVF is Intel Visual Fortran, for which the IMSL library is planned but not yet available. === Subject: Re: Examples using IMSL for non-linear curve fitting The IMSL subroutines UNLSF and UNLSJ can be used to solve a nonlinear least-squares problem using a modified Levenberg-Marquardt algorithm. See Shapter 8 of their Math/Library manual for descriptions and examples of calling these subroutines. > > The IMSL for FORTRAN library has routines to perform non-linear least square > curve fitting. Is there any example FORTRAN code using IMSL for non-linear > curve fitting available so I can learn how to use the library. I will use > CVF 6.6. > > > The user guides at > include examples of use. === Subject: Re: Round off example: MATLAB ok, GCC not ok I have discussed this on a Norwegian group and I know now what the problem is: It depends on how many bits that are used for the floating point registers. As you know, the 80 bits. But the OS can choose to only use 64 of these. Linux uses all 80 bits and gets stable result, FreeBSD uses 64 bits (IEEE) and gets unstable results. Matlab uses 64 bits and gets instable results. Windows uses 64 bits and gets instabile results (with MSVC or Microsoft's gcc build, which uses Windows' standard behaviour by default.) However, MinGW overrides Windows standard behaviour and turns the resolution up to 80 bits like Linux, and the calculation then becomes stable. (Cygwin behaves like MinGW.) This explains all the variability between environments systems and compilers. If you use GCC on Linux or MinGW on Windows, you can explicitely set the floating point resolution down to 64 bits using the instruction _controlfp(_PC_64, _MCW_PC); which is declared in header file . This forces unstable behaviour even on Linux or Windows with MinGW. You can also get unstable result by overflowing the 80 bit resolution, e.g. by using z = 4000.0*x - 3999.0*y; Sturla Molden === Subject: Re: Round off example: MATLAB ok, GCC not ok > The algorithm below is unstable under > MATLAB since the number 2.1 > has an infinite binary expansion. ... > format long > > x = 2.1; > y = 2.1; > > for i = 1:40 = (110).(1001)^{10}(1001)( > z = 4.0*x - 3.0*y > x = y; > y = z; > end Are you sure? That's not how I figure it. Let's try writing everything in binary and using IEEE double-precision arithmetic (so, 52-bit mantissas, or 52 digits after the initial 1). 2.1 = (in binary) 10.(0011) recurring represented in IEEE double-precision arithmetic as (10).(0011)^{12}(010) which is the initial value of x & y (^{12] just means repeated (decimal) 12 times). Then in the first loop 4x = (1000).(1100)^{11}(1101) 3y = (110).(1001)^{12}(110) which gets rounded to (110).(1001)^{12}(11) = (110).(1001)^{111}(1001)(11) Subtracting the two gives us 4x - 3y = (10).(0011)^{11}(0011)(01) (we have to borrow one) = (10).(0011)^{12}(01) meaning z is the same as the initial value of x & y. So it looks to me as if GCC is right and Matlab (assuming it claims to implement IEEE double precision arithmetic throughout this calculation) is wrong. === Subject: Re: Round off example: MATLAB ok, GCC not ok No sorry, that's balderdash! Please ignore!! === Subject: Re: Round off example: MATLAB ok, GCC not ok No hope! I've tried to use the ``volatile'' directive, to put ``4.0*x-3.0*y'' inside a function as well to turn off compiler optimization. I still get the same result. I even tried to compile the C code using the old Turbo C 2.01 compiler. No way. :( Humberto. > > The algorithm below is unstable under > MATLAB since the number 2.1 > has an infinite binary expansion. > > I've tried to reproduce the same behaviour > under GCC but, no luck, it's stable with > the famous C compiler. > > What's happening? Is it possible > to make GCC to act like MATLAB? Am I missing > some compiler directive? > > > --------------------------------------- > format long > > x = 2.1; > y = 2.1; > > for i = 1:40 > z = 4.0*x - 3.0*y > x = y; > y = z; > end > --------------------------------------- > > The Gnu compiler is pretty agressive about optimizing -- I could believe > that it would see that x and y are equal, and therefore that z = y already, > and basically optimize your whole loop to nothing. > > You can try to turn optimizations off, but I've had that compiler optimize > some anyway (maybe it's streamlining, not optimization). > > If you want to absolutely _force_ it to make the computation put the 4.0 * > x - 3.0 * y into a function call, or declare x and y to be volatile (the > function call would be a stronger way to do it). > > > ---------------------------------------- > Tim Wescott > Wescott Design Services > http://www.wescottdesign.com === Subject: Re: Round off example: MATLAB ok, GCC not ok > The Gnu compiler is pretty agressive about optimizing -- I could believe > that it would see that x and y are equal, and therefore that z = y already, > and basically optimize your whole loop to nothing. I'm note sure what is happening. gcc is behaving strangely. Let us try the C code: #include int main() { double x = 2.1, y = 2.1, z; int i; for (i=0; i<40; i++) { z = 4.0*x - 3.0*y; x = y; y = z; } printf(x = %f,ty = %f,n,x,y); } First I test with the MinGW distribution of GCC 3.2 on Windows XP, using the MSYS bash shell: $ gcc -o test.exe main.c $ ./test x = 2.100000, y = 2.100000, Then turning the optimization off does not help: $ gcc -O0 -o test.exe main.c $ ./test x = 2.100000, y = 2.100000, Declaring x,y and z as volatile does not help $ gcc -O0 -o test.exe main.c $ ./test x = 2.100000, y = 2.100000, Setting the loop counter i to volatile does not help: $ gcc -O0 -o test.exe main.c $ ./test x = 2.100000, y = 2.100000, Setting x, y, z, and i to volatile does not help: $ gcc -O0 -o test.exe main.c $ ./test x = 2.100000, y = 2.100000, However, using Microsoft's build of GCC 3.3 (the gcc that ships with Windows Services for Unix 3.5) does help. With no volatile variables and different levels of optimization: % gcc -o test main.c % ./test x = -26843543.500000, y = 107374184.500000, % gcc -O2 -o test main.c % ./test x = -26843543.500000, y = 107374184.500000, % gcc -O3 -o test main.c % ./test x = -26843543.500000, y = 107374184.500000, % Microsoft's own compiler also works. With full optimization and no volatile variables: C:developmentgnutest>cl -O2 -o test.exe main.c Microsoft (R) 32-bit C/C++ Optimizing Compiler Version 13.10.3077 for 80x86 Copyright (C) Microsoft Corporation 1984-2002. All rights reserved. main.c Microsoft (R) Incremental Linker Version 7.10.3077 Copyright (C) Microsoft Corporation. All rights reserved. /out:main.exe /out:test.exe main.obj C:developmentgnutest>test x = -26843543.500000, y = 107374184.500000, Then finally Matlab, just to prove that it works: function gnutest x = 2.1; y = 2.1; for i = 1:40 z = 4.0*x - 3.0*y; x = y; y = z; end fprintf(1,'x = %f,ty = %f,tz = %fn',z,y,z); >> gnutest x = 107374184.500000, y = 107374184.500000, z = 107374184.500000 As it should be. Sturla Molden === Subject: Re: Round off example: MATLAB ok, GCC not ok >#include >int main() >{ > double x = 2.1, y = 2.1, z; > int i; > for (i=0; i<40; i++) > { > z = 4.0*x - 3.0*y; > x = y; > y = z; > } > printf(x = %f,ty = %f,n,x,y); >} The behavior of this code with different compilers on Intel architecture machines can be explained by whether or not the 80x87's extended precision registers are being used. I compiled the code with Intel's C compiler, first using -pc80 (uses extended precision registers when possible) and then with -pc64 (all results will be rounded to double precision), and got: >icc -pc80 test.c -o test80 >icc -pc64 test.c -o test64 >./test64 x = -26843543.500000, y = 107374184.500000, >./test80 x = 2.100000, y = 2.100000, Every version of gcc that I have used under linux uses the extended precision registers unless you tell it to use SSE2 instructions (the SSE2 instructions use a different set of floating point registers that only hold single or double precision numbers.) You can force gcc to produce code that will round all results to double precision by using the -ffloat-store compiler switch. A better alternative if you want to force double precision instead of extended precision is to set the precision control to 64 bits at the start of the program. Note that if you ran this C code on a machine that didn't have extended precision (say a Sun Ultrasparc machine, or an Apple Macintosh with a G4 or G5 processor), then you would normally get results like the output from test64. Note also that MATLAB always runs in 64 bit precision, even on machines with support for extended precision. -- Brian Borchers borchers@nmt.edu Department of Mathematics http://www.nmt.edu/~borchers/ New Mexico Tech Phone: 505-835-5813 Socorro, NM 87801 FAX: 505-835-5366 http://www.newsfeeds.com - The #1 Newsgroup Service in the World! -----== Over 100,000 Newsgroups - 19 Different Servers! =----- === Subject: Re: Round off example: MATLAB ok, GCC not ok >> The Gnu compiler is pretty agressive about optimizing -- I could believe >> that it would see that x and y are equal, and therefore that z = y already, >> and basically optimize your whole loop to nothing. > I'm note sure what is happening. gcc is behaving strangely. > Let us try the C code: > #include > int main() > { > double x = 2.1, y = 2.1, z; > int i; > for (i=0; i<40; i++) > { > z = 4.0*x - 3.0*y; > x = y; > y = z; > } > printf(x = %f,ty = %f,n,x,y); > } I've been staring at this for while, and I think I know what's happening: `-ffloat-store' Do not store floating point variables in registeand inhibit other options that might change whether a floating point value is taken from a register or memory. This option prevents undesirable excess precision on machines such as the 68000 where the floating registers (of the 68881) keep more precision than a `double' is supposed to have. Similarly for the x86 architecture. For most programs, the excess precision does only good, but a few programs rely on the precise definition of IEEE floating point. Use `-ffloat-store' for such programs, after modifying them to store all pertinent intermediate computations into variables. The Intel architecure has 80-bit floating point registeinstead of the 64-bit specified for IEEE doubles. As the GCC text says, that's usually good. I was able to reproduce the unstable results by computing intermediate results for 4.0*x and 3.0*y, then subtracting. The optimizer will kill those intermediate computations unless you specify -ffloat-store. The other choice is to use the SSE unit (on Pentium 4 processors), which is faster for those computations it supports, and has only 64 bit registers. Using GCC, specify these options to use SSE math, and original program will be unstable: -march=pentium4 -mfpmath=sse So I think this all has to do with which compilers/environments a) Store intermediate results in variables b) Use SSE when supported MATLAB falls into a), since the matlab function calls return values in memory. Sturla's other examples of Microsoft's GCC build and MSVC probably have SSE math enabled by default, for speed. -- Peter Boettcher MIT Lincoln Laboratory MATLAB FAQ: http://www.mit.edu/~pwb/cssm/ === Subject: Re: Round off example: MATLAB ok, GCC not ok Right. MATLAB mimics 64 bit floating point registers by setting the precision of the chip to a 53 bit mantissa (including the implicit leading 1) even when A c compiler gets to keep a simple computation like z = 4.0*x - 3.0*y; in the 80 bit register as long as it can but may store some of it out to 64 bit double memory depending. You can see exactly the choices made by the compiler in the generated assembly code. If MATLAB used the full extended precision of these chips (64 bit mantissa) in combination with the JIT optimizations that no longer force the code snippet to become: tmp1 = 4.0 * x; % 64 bit double temporary storage tmp2 = 3.0 * y; % bit bit double temporary storage z = tmp1 - tmp2; % 64 bit double results you would see the same results as gcc. Penny. -- Penny Anderson The MathWorks Please post any followups to this newsgroup. > >> The Gnu compiler is pretty agressive about optimizing -- I could believe >> that it would see that x and y are equal, and therefore that z = y already, >> and basically optimize your whole loop to nothing. > > I'm note sure what is happening. gcc is behaving strangely. > Let us try the C code: > > #include > int main() > { > double x = 2.1, y = 2.1, z; > int i; > for (i=0; i<40; i++) > { > z = 4.0*x - 3.0*y; > x = y; > y = z; > } > printf(x = %f,ty = %f,n,x,y); > } > I've been staring at this for while, and I think I know what's > happening: > `-ffloat-store' > Do not store floating point variables in registeand inhibit > other options that might change whether a floating point value is > taken from a register or memory. > This option prevents undesirable excess precision on machines such > as the 68000 where the floating registers (of the 68881) keep more > precision than a `double' is supposed to have. Similarly for the > x86 architecture. For most programs, the excess precision does > only good, but a few programs rely on the precise definition of > IEEE floating point. Use `-ffloat-store' for such programs, after > modifying them to store all pertinent intermediate computations > into variables. > The Intel architecure has 80-bit floating point registeinstead of > the 64-bit specified for IEEE doubles. As the GCC text says, that's > usually good. I was able to reproduce the unstable results by > computing intermediate results for 4.0*x and 3.0*y, then subtracting. > The optimizer will kill those intermediate computations unless you > specify -ffloat-store. > The other choice is to use the SSE unit (on Pentium 4 processors), > which is faster for those computations it supports, and has only 64 > bit registers. Using GCC, specify these options to use SSE math, and > original program will be unstable: -march=pentium4 -mfpmath=sse > So I think this all has to do with which compilers/environments > a) Store intermediate results in variables > b) Use SSE when supported > MATLAB falls into a), since the matlab function calls return values in > memory. > Sturla's other examples of Microsoft's GCC build and MSVC probably > have SSE math enabled by default, for speed. > -- > Peter Boettcher > MIT Lincoln Laboratory > MATLAB FAQ: http://www.mit.edu/~pwb/cssm/ === Subject: Re: Round off example: MATLAB ok, GCC not ok > MATLAB mimics 64 bit floating point registers by setting the > precision of the chip to a 53 bit mantissa (including the > implicit leading 1) even when 80 bit registers are used on Which, by the way, is the topic of an interesting paper by William Kahan: http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf Peter -- If I made a copy of the The Digital Millennium Copyright Act, would that be a violation of it? === Subject: Re: Round off example: MATLAB ok, GCC not ok > Sturla's other examples of Microsoft's GCC build and MSVC probably > have SSE math enabled by default, for speed. I'm not sure. I asked on a Norwegian newsgroup and one person reported different results with the same executable on the same CPU. The code was compiled with gcc on Linux and then the binary was tested on Linux and FreeBSD (using linux emulation). On Linux the computation was reported stabile, on FreeBSD the computation was reported unstabile. Inside the Intel processors the SSE and floating point registers are separate. The MMX registers are 64 bits and occupy the same space as the 80 bit floating point registers; the SSE registers are 128 bit and is a separate entity in the CPU. Floating point operations is compiled using the floating point registers. Non-standard code e.g. compiler pragmas, inline assembly or streaming SIMD intrinsics are required to touch the SSE registers. I checked the assembly code produced by the two gcc builds I tested (by compiling with gcc -S). The assembly code was essentially identical, and no SSE opcode was used. Still, the Win32 and Interix subsytems on my computer produced different results, just like Linux and FreeBSD with the same binary. Sturla Molden BTW: I would like to check the assembly from MSVC as well, but I cannot find a compiler switch to output the assembly code. === Subject: Finite Volume Implementation! by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i29HbmL15269; Hi there! implementation of the finite volume method is fully explained? Or where can I find an example program with the appropriate documentation? George === Subject: How to solve this differention equation? I encountered a differention equation problem as below. Assume g(x) is a a differentiable scalar function of x, and the equation is (g'(x))^2 = 1/(1+g^2(x)), where g'(x) is the first order direvative of g(x) w.r.t x. So given the above equation, how to get the explicit form of g(x) or implicit value of g(x) give a particular value x. Fred === Subject: Re: How to solve this differention equation? > >I encountered a differention equation problem as below. > >Assume g(x) is a a differentiable scalar function of x, >and the equation is >(g'(x))^2 = 1/(1+g^2(x)), where g'(x) is the first order direvative of >g(x) w.r.t x. > >So given the above equation, how to get the >explicit form of g(x) or implicit value of g(x) give a particular value x. > you have two possible equations: g'(x) = -1/sqrt(1+g(x)^2) or g'(x) = 1/sqrt(1+g(x)^2) using separation of variables you get pm integral{from y0 to y} sqrt(1+u^2)du = x-x0 were pm is + or - pm ( (y/2)sqrt(1+y^2)+(1/2)*log( y+sqrt(1+y^2)) - (y0/2)sqrt(1+y0^2) - (1/2)*log(y0+sqrt(1+y0^2) )) = x-x0; to be solved for y=y(x) given the initial value (x0,y0) hth peter === Subject: demande by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i29HtrW17714; Salut j'ai l'honneur de vous demand.8e de m'envoiyer un programme qui fait l'integration en math.8ematique. Dans l'attente d'une r.8eponse favorable, veiller acc.8epter mes expressions les plus distingu.8ees. === Subject: Runge-Kutta Fourth Order Method in Maple code by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i29Kke507322; I have to solve van der pols equation by the runge-kutta fourth order method, by programming the method in maple, However I have been told that maple has the fourth order method built in and to program it in maple is trivial. I am awful at maple and really need to know if any one has the code they could e-mail me PLEASE!! Also, it would be very helpful if i could get the results in a table and be able to plot them. PLEASE -I REALLY NEED THE CODE. THANKS IN ADVANCE === Subject: Uw vraag by support1.mathforum.org (8.11.6/8.11.6/The Math Forum, $Revision: 1.9 primary) id i29MJJV18946; Ik weet het ook niet.