rewrote rotmg.c instead of modifying very old code

Author: wernsaar

interface/rotmg.c

@@ -1,3 +1,37 @@
+/***************************************************************************
+Copyright (c) 2013, The OpenBLAS Project
+All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+1. Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in
+the documentation and/or other materials provided with the
+distribution.
+3. Neither the name of the OpenBLAS project nor the names of
+its contributors may be used to endorse or promote products
+derived from this software without specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
+USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*****************************************************************************/
+
+/**************************************************************************************
+* 2014/02/28 Saar
+*	Test with lapack-3.5.0	: OK 
+*
+**************************************************************************************/
+
+
 #include "common.h"
 #ifdef FUNCTION_PROFILE
 #include "functable.h"
@@ -7,6 +41,8 @@
 #define  GAMSQ   16777216.e0
 #define  RGAMSQ  5.9604645e-8
 
+#define  TWO 2.e0
+
 #ifdef DOUBLE
 #define ABS(x) fabs(x)
 #else
@@ -25,181 +61,168 @@ void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
 
 #endif
 
-    FLOAT du, dp1, dp2, dq2, dq1, dh11, dh21, dh12, dh22;
-    int igo, flag;
-    FLOAT dtemp;
-
-#ifndef CBLAS
-  PRINT_DEBUG_NAME;
-#else
-  PRINT_DEBUG_CNAME;
-#endif
-
-    dh11 = ZERO;
-    dh12 = ZERO;
-    dh21 = ZERO;
-    dh22 = ZERO;
-
-    if (*dd1 < ZERO) goto L60;
-
-    dp2 = *dd2 * dy1;
-
-    if (dp2 == ZERO) {
-      flag = -2;
-      goto L260;
-    }
-
-    dp1 = *dd1 * *dx1;
-    dq2 =  dp2 * dy1;
-    dq1 =  dp1 * *dx1;
-
-    if (! (ABS(dq1) > ABS(dq2))) goto L40;
-
-    dh21 = -(dy1) / *dx1;
-    dh12 = dp2 / dp1;
-
-    du = ONE - dh12 * dh21;
-
-    if (du <= ZERO) goto L60;
-
-    flag = 0;
-    *dd1 /= du;
-    *dd2 /= du;
-    *dx1 *= du;
-
-    goto L100;
-
-L40:
-    if (dq2 < ZERO) goto L60;
-
-    flag = 1;
-    dh11  = dp1 / dp2;
-    dh22  = *dx1 / dy1;
-    du    = ONE + dh11 * dh22;
-    dtemp = *dd2 / du;
-    *dd2  = *dd1 / du;
-    *dd1  = dtemp;
-    *dx1  = dy1 * du;
-    goto L100;
-
-L60:
-    flag = -1;
-    dh11 = ZERO;
-    dh12 = ZERO;
-    dh21 = ZERO;
-    dh22 = ZERO;
-
-    *dd1 = ZERO;
-    *dd2 = ZERO;
-    *dx1 = ZERO;
-    goto L220;
-
-
-L70:
-    if (flag < 0) goto L90;
- 
-    if (flag > 0) goto L80;
- 
-    dh11 = ONE;
-    dh22 = ONE;
-    flag = -1;
-    goto L90;
-
-L80:
-    dh21 = -ONE;
-    dh12 = ONE;
-    flag = -1;
-
-L90:
-    switch (igo) {
-	case 0: goto L120;
-	case 1: goto L150;
-	case 2: goto L180;
-	case 3: goto L210;
-    }
-
-L100:
-    if (!(*dd1 <= RGAMSQ)) goto L130;
-    if (*dd1 == ZERO) goto L160;
-    igo = 0;
-    goto L70;
-
-L120:
-    *dd1 *= GAM * GAM;
-    *dx1 /= GAM;
-    dh11 /= GAM;
-    dh12 /= GAM;
-    goto L100;
-
-L130:
-    if (! (*dd1 >= GAMSQ)) {
-	goto L160;
-    }
-    igo = 1;
-    goto L70;
-
-L150:
-    *dd1 /= GAM * GAM;
-    *dx1 *= GAM;
-    dh11 *= GAM;
-    dh12 *= GAM;
-    goto L130;
-
-L160:
-    if (! (ABS(*dd2) <= RGAMSQ)) {
-	goto L190;
-    }
-    if (*dd2 == ZERO) {
-	goto L220;
-    }
-    igo = 2;
-    goto L70;
-
-L180:
-/* Computing 2nd power */
-    *dd2 *= GAM * GAM;
-    dh21 /= GAM;
-    dh22 /= GAM;
-    goto L160;
-
-L190:
-    if (! (ABS(*dd2) >= GAMSQ)) {
-	goto L220;
-    }
-    igo = 3;
-    goto L70;
-
-L210:
-/* Computing 2nd power */
-    *dd2 /= GAM * GAM;
-    dh21 *= GAM;
-    dh22 *= GAM;
-    goto L190;
-
-L220:
-    if (flag < 0) {
-	goto L250;
-    } else if (flag == 0) {
-	goto L230;
-    } else {
-	goto L240;
-    }
-L230:
-    dparam[2] = dh21;
-    dparam[3] = dh12;
-    goto L260;
-L240:
-    dparam[2] = dh11;
-    dparam[4] = dh22;
-    goto L260;
-L250:
-    dparam[1] = dh11;
-    dparam[2] = dh21;
-    dparam[3] = dh12;
-    dparam[4] = dh22;
-L260:
-    dparam[0] = (FLOAT) flag;
-    return;
+	FLOAT du, dp1, dp2, dq2, dq1, dh11, dh21, dh12, dh22, dflag, dtemp;
+
+	if(*dd1 < ZERO)
+	{
+		dflag = -ONE;
+		dh11  = ZERO;
+		dh12  = ZERO;
+		dh21  = ZERO;
+		dh22  = ZERO;
+
+		*dd1  = ZERO;
+		*dd2  = ZERO;
+		*dx1  = ZERO;
+	}
+	else
+	{
+		dp2 = *dd2 * dy1;
+		if(dp2 == ZERO)
+		{
+			dflag = -TWO;
+			dparam[0] = dflag;
+			return;
+		}
+		dp1 = *dd1 * *dx1;
+		dq2 =  dp2 * dy1;
+		dq1 =  dp1 * *dx1;
+		if(ABS(dq1) > ABS(dq2))
+		{
+			dh21 = -  dy1 / *dx1;
+			dh12 =    dp2 /  dp1;
+
+			du   = ONE - dh12 * dh21;
+			if(du > ZERO)
+			{
+				dflag = ZERO;
+				*dd1  = *dd1 / du;
+				*dd2  = *dd2 / du;
+				*dx1  = *dx1 * du;
+
+			}
+		}
+		else
+		{
+			if(dq2 < ZERO)
+			{
+				dflag = -ONE;
+				
+				dh11  = ZERO;
+				dh12  = ZERO;
+				dh21  = ZERO;
+				dh22  = ZERO;
+
+				*dd1  = ZERO;
+				*dd2  = ZERO;
+				*dx1  = ZERO;
+			}
+			else
+			{
+				dflag = ONE;
+
+				dh11  =  dp1 /  dp2;
+				dh22  = *dx1 /  dy1;
+				du    =  ONE + dh11 * dh22;
+				dtemp = *dd2 / du;
+
+				*dd2  = *dd1 / du;
+				*dd1  = dtemp;
+				*dx1  = dy1 * du;
+			}			
+		}
+
+
+		if(*dd1 != ZERO)
+		{
+			while( (*dd1 <= RGAMSQ) || (*dd1 >= GAMSQ) )
+			{
+				if(dflag == ZERO)
+				{
+					dh11  =  ONE;
+					dh22  =  ONE;
+					dflag = -ONE;
+				}
+				else
+				{
+					dh21  = -ONE;
+					dh12  =  ONE;
+					dflag = -ONE;
+				}
+				if( *dd1 <= RGAMSQ )
+				{
+					*dd1  = *dd1 * (GAM * GAM);
+					*dx1  = *dx1 / GAM;
+					dh11  = dh11 / GAM;
+					dh12  = dh12 / GAM;
+				}
+				else
+				{
+					*dd1  = *dd1 / (GAM * GAM);
+					*dx1  = *dx1 * GAM;
+					dh11  = dh11 * GAM;
+					dh12  = dh12 * GAM;
+				}
+			}
+		}
+			
+		if(*dd2 != ZERO)
+		{
+			while( (ABS(*dd2) <= RGAMSQ) || (ABS(*dd2) >= GAMSQ) )
+			{
+				if(dflag == ZERO)
+				{
+					dh11  =  ONE;
+					dh22  =  ONE;
+					dflag = -ONE;
+				}
+				else
+				{
+					dh21  = -ONE;
+					dh12  =  ONE;
+					dflag = -ONE;
+				}
+				if( ABS(*dd2) <= RGAMSQ )
+				{
+					*dd2  = *dd2 * (GAM * GAM);
+					dh21  = dh21 / GAM;
+					dh22  = dh22 / GAM;
+				}
+				else
+				{
+					*dd2  = *dd2 / (GAM * GAM);
+					dh21  = dh21 * GAM;
+					dh22  = dh22 * GAM;
+				}
+			}
+		}
+					
+	}
+
+	if(dflag < ZERO)
+	{
+		dparam[1] = dh11;
+		dparam[2] = dh21;
+		dparam[3] = dh12;
+		dparam[4] = dh22;
+	}
+	else
+	{
+		if(dflag == ZERO)
+		{
+			dparam[2] = dh21;
+			dparam[3] = dh12;
+		}
+		else
+		{
+			dparam[1] = dh11;
+			dparam[4] = dh22;
+		}
+	}
+
+
+	dparam[0] = dflag;
+	return;
 }