interface/rotmg.c

   1 /***************************************************************************
   2 Copyright (c) 2013, The OpenBLAS Project
   3 All rights reserved.
   4 Redistribution and use in source and binary forms, with or without
   5 modification, are permitted provided that the following conditions are
   6 met:
   7 1. Redistributions of source code must retain the above copyright
   8 notice, this list of conditions and the following disclaimer.
   9 2. Redistributions in binary form must reproduce the above copyright
  10 notice, this list of conditions and the following disclaimer in
  11 the documentation and/or other materials provided with the
  12 distribution.
  13 3. Neither the name of the OpenBLAS project nor the names of
  14 its contributors may be used to endorse or promote products
  15 derived from this software without specific prior written permission.
  16 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  17 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  18 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  19 ARE DISCLAIMED. IN NO EVENT SHALL THE OPENBLAS PROJECT OR CONTRIBUTORS BE
  20 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  21 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  22 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  23 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  24 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
  25 USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  26 *****************************************************************************/
  27
  28 /**************************************************************************************
  29 * 2014/05/02 Saar
  30 *       fixed two bugs as reported by Brendan Tracey
  31 *       Test with lapack-3.5.0  : OK
  32 *
  33 **************************************************************************************/
  34
  35
  36 #include "common.h"
  37 #ifdef FUNCTION_PROFILE
  38 #include "functable.h"
  39 #endif
  40
  41 #define  GAM     4096.e0
  42 #define  GAMSQ   16777216.e0
  43 #define  RGAMSQ  5.9604645e-8
  44
  45 #define  TWO 2.e0
  46
  47 #ifdef DOUBLE
  48 #define ABS(x) fabs(x)
  49 #else
  50 #define ABS(x) fabsf(x)
  51 #endif
  52
  53 #ifndef CBLAS
  54
  55 void NAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT *DY1, FLOAT *dparam){
  56
  57   FLOAT dy1 = *DY1;
  58
  59 #else
  60
  61 void CNAME(FLOAT *dd1, FLOAT *dd2, FLOAT *dx1, FLOAT dy1, FLOAT *dparam){
  62
  63 #endif
  64
  65         FLOAT du, dp1, dp2, dq2, dq1, dh11=ZERO, dh21=ZERO, dh12=ZERO, dh22=ZERO, dflag=-ONE, dtemp;
  66
  67         if (*dd2 == ZERO || dy1 == ZERO)
  68         {
  69                 dflag = -TWO;
  70                 dparam[0] = dflag;
  71                 return;
  72         }
  73
  74         if(*dd1 < ZERO)
  75         {
  76                 dflag = -ONE;
  77                 dh11  = ZERO;
  78                 dh12  = ZERO;
  79                 dh21  = ZERO;
  80                 dh22  = ZERO;
  81
  82                 *dd1  = ZERO;
  83                 *dd2  = ZERO;
  84                 *dx1  = ZERO;
  85         }
  86         else if ((*dd1 == ZERO || *dx1 == ZERO) && *dd2 > ZERO)
  87         {
  88                 dflag = ONE;
  89                 dh12 = 1;
  90                 dh21 = -1;
  91                 *dx1 = dy1;
  92                 dtemp = *dd1;
  93                 *dd1 = *dd2;
  94                 *dd2 = dtemp;
  95         }
  96         else
  97         {
  98                 dp2 = *dd2 * dy1;
  99                 if(dp2 == ZERO)
 100                 {
 101                         dflag = -TWO;
 102                         dparam[0] = dflag;
 103                         return;
 104                 }
 105                 dp1 = *dd1 * *dx1;
 106                 dq2 =  dp2 * dy1;
 107                 dq1 =  dp1 * *dx1;
 108                 if(ABS(dq1) > ABS(dq2))
 109                 {
 110                         dflag = ZERO;
 111                         dh11  =  ONE;
 112                         dh22  =  ONE;
 113                         dh21 = -  dy1 / *dx1;
 114                         dh12 =    dp2 /  dp1;
 115
 116                         du   = ONE - dh12 * dh21;
 117                         if(du > ZERO)
 118                         {
 119                                 dflag = ZERO;
 120                                 *dd1  = *dd1 / du;
 121                                 *dd2  = *dd2 / du;
 122                                 *dx1  = *dx1 * du;
 123                         } else {
 124                                 dflag = -ONE;
 125
 126                                 dh11  = ZERO;
 127                                 dh12  = ZERO;
 128                                 dh21  = ZERO;
 129                                 dh22  = ZERO;
 130
 131                                 *dd1  = ZERO;
 132                                 *dd2  = ZERO;
 133                                 *dx1  = ZERO;
 134                         }
 135
 136                 }
 137                 else
 138                 {
 139                         if(dq2 < ZERO)
 140                         {
 141                                 dflag = -ONE;
 142
 143                                 dh11  = ZERO;
 144                                 dh12  = ZERO;
 145                                 dh21  = ZERO;
 146                                 dh22  = ZERO;
 147
 148                                 *dd1  = ZERO;
 149                                 *dd2  = ZERO;
 150                                 *dx1  = ZERO;
 151                         }
 152                         else
 153                         {
 154                                 dflag =  ONE;
 155                                 dh21  = -ONE;
 156                                 dh12  =  ONE;
 157
 158                                 dh11  =  dp1 /  dp2;
 159                                 dh22  = *dx1 /  dy1;
 160                                 du    =  ONE + dh11 * dh22;
 161                                 dtemp = *dd2 / du;
 162
 163                                 *dd2  = *dd1 / du;
 164                                 *dd1  = dtemp;
 165                                 *dx1  = dy1 * du;
 166                         }
 167                 }
 168
 169
 170                 while ( *dd1 <= RGAMSQ && *dd1 != ZERO)
 171                 {
 172                         dflag = -ONE;
 173                         *dd1  = *dd1 * (GAM * GAM);
 174                         *dx1  = *dx1 / GAM;
 175                         dh11  = dh11 / GAM;
 176                         dh12  = dh12 / GAM;
 177                 }
 178                 while (ABS(*dd1) > GAMSQ) {
 179                         dflag = -ONE;
 180                         *dd1  = *dd1 / (GAM * GAM);
 181                         *dx1  = *dx1 * GAM;
 182                         dh11  = dh11 * GAM;
 183                         dh12  = dh12 * GAM;
 184                 }
 185
 186                 while (ABS(*dd2) <= RGAMSQ && *dd2 != ZERO) {
 187                         dflag = -ONE;
 188                         *dd2  = *dd2 * (GAM * GAM);
 189                         dh21  = dh21 / GAM;
 190                         dh22  = dh22 / GAM;
 191                 }
 192                 while (ABS(*dd2) > GAMSQ) {
 193                         dflag = -ONE;
 194                         *dd2  = *dd2 / (GAM * GAM);
 195                         dh21  = dh21 * GAM;
 196                         dh22  = dh22 * GAM;
 197                 }
 198
 199         }
 200
 201         if(dflag < ZERO)
 202         {
 203                 dparam[1] = dh11;
 204                 dparam[2] = dh21;
 205                 dparam[3] = dh12;
 206                 dparam[4] = dh22;
 207         }
 208         else
 209         {
 210                 if(dflag == ZERO)
 211                 {
 212                         dparam[2] = dh21;
 213                         dparam[3] = dh12;
 214                 }
 215                 else
 216                 {
 217                         dparam[1] = dh11;
 218                         dparam[4] = dh22;
 219                 }
 220         }
 221
 222
 223         dparam[0] = dflag;
 224         return;
 225 }
 226
 227