glyman
04-29-2005, 10:27 PM
I recently had Powell cause an unpleasant crash. While the book discusses the problem of what to do when no progress is made in one of the directions (xit has zero norm), I recently added a fix to avoid the problem.
Powell's normal method of getting a new set of directions (columns of array 'xi') is to replace column 'ibig' with the last or nth column and the displacement vector 'xit' is then written into the last column of 'xi'. This is great as long as 'xit' is not a null vector (which can happen!!). If it is null, the direction set 'xi' no longer spans the parameter space and mnbrak can be asked to search in a null direction. this is cause and overflow in mnbrak as it seaches to 'infinity' looking for a better point.
The fix is to check the norm of 'xit' (sum of squares of the elements will do) and if it is zero, do an svd on the xi matrix that has had the null vector inserted. You can then take the 'U' matrix that comes out of the svd and rescxale it columns in line with the columns of xi before the svd.
The following bit of code inserted into the 'do 15' loop (last do loop in hte Powell routine will do the trick.
norm = 0.0d0
do i=1,n
norm = norm*xit(i)**2
enddo
if(norm .eq. 0.0d0) then
allocate(xfix(n,n),vfix(n,n),wfix(n))
xfix = xi
call svdcmp(xfix,n,n,n,n,wfix,vfix)
wfix(n) = minval(wfix(1:n-1))
do i=1,n-1
xi(:,i) = xfix(:,i)*dot_product(xi(:,i),xfix(:,i))
enddo
xi(:,n) = wfix(n)*xfix(:,n)
endif
don't forget to deallocate the arrays!!
Powell's normal method of getting a new set of directions (columns of array 'xi') is to replace column 'ibig' with the last or nth column and the displacement vector 'xit' is then written into the last column of 'xi'. This is great as long as 'xit' is not a null vector (which can happen!!). If it is null, the direction set 'xi' no longer spans the parameter space and mnbrak can be asked to search in a null direction. this is cause and overflow in mnbrak as it seaches to 'infinity' looking for a better point.
The fix is to check the norm of 'xit' (sum of squares of the elements will do) and if it is zero, do an svd on the xi matrix that has had the null vector inserted. You can then take the 'U' matrix that comes out of the svd and rescxale it columns in line with the columns of xi before the svd.
The following bit of code inserted into the 'do 15' loop (last do loop in hte Powell routine will do the trick.
norm = 0.0d0
do i=1,n
norm = norm*xit(i)**2
enddo
if(norm .eq. 0.0d0) then
allocate(xfix(n,n),vfix(n,n),wfix(n))
xfix = xi
call svdcmp(xfix,n,n,n,n,wfix,vfix)
wfix(n) = minval(wfix(1:n-1))
do i=1,n-1
xi(:,i) = xfix(:,i)*dot_product(xi(:,i),xfix(:,i))
enddo
xi(:,n) = wfix(n)*xfix(:,n)
endif
don't forget to deallocate the arrays!!