The Perils of Polishing (LONG!)

            do while (cabs1(det(1)) .ge. ten)
               det(1) = det(1) * 0.1e0
               det(2) = det(2) + 1.0e0
            end do

I have seen code like this (or the f77 version) in several linear algebra libraries. The problem is that determinants can span many orders of magnitude and can easily exceed the floating point exponent range. These loops keep track separately of the exponent and the scaled value, avoiding the overflow and underflow problem. The posted code is for complex determinants, but the same issue arises for real determinants.

Another approach to handle this situation involves computing the log of the determinant rather than the determinant itself. The corner case is when the determinant is zero, which must be handled separately if it can occur. If the sign or complex phase is arbitrary, or if the intermediate signs/phases are important, then that must be tracked separately too.

In modern fortran, I think I would rather do this computation with a derived type. One member would be the det(1) value, i.e. a scaled value that is kept within some restricted range of values. The other member would be the integer exponent, stored above in det(2). Of course, if working with a legacy code, then this argument change would need to be propagated throughout.

I’ve noticed that the logic is almost always done with decimal exponents. This has always puzzled me. I would think that it would be more natural to use binary exponent ranges. That way, all of the exponent scaling can be done by just shifting the exponent values, there would be no need to actually do the floating point multiplications and divisions or to touch the mantissa. In modern fortran, there are the EXPONENT() and SCALE() intrinsics which allow this to be done in a portable way.