OpenBLAS/interface/zrotg.c

#include <math.h>
#include <float.h>
#include "common.h"
#ifdef FUNCTION_PROFILE
#include "functable.h"
#endif


#ifndef CBLAS
void NAME(FLOAT *DA, FLOAT *DB, FLOAT *C, FLOAT *S){

#else
void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {
    FLOAT *DA = (FLOAT*) VDA;
    FLOAT *DB = (FLOAT*) VDB;
    FLOAT *S  = (FLOAT*) VS;
#endif /* CBLAS */

#ifdef DOUBLE
  long double safmin = DBL_MIN;
#else
  long double safmin = FLT_MIN;
#endif

#if defined(__i386__) || defined(__x86_64__) || defined(__ia64__) || defined(_M_X64) || defined(_M_IX86)

  long double da_r = *(DA + 0);
  long double da_i = *(DA + 1);
  long double db_r = *(DB + 0);
  long double db_i = *(DB + 1);
  long double r;

  long double ada = fabsl(da_r) + fabsl(da_i);
  long double adb = sqrt(db_r * db_r + db_i * db_i);

  PRINT_DEBUG_NAME;

  IDEBUG_START;

  FUNCTION_PROFILE_START();

  if (ada == ZERO) {
    *C        = ZERO;
    *(S  + 0) = ONE;
    *(S  + 1) = ZERO;
    *(DA + 0) = db_r;
    *(DA + 1) = db_i;
  } else {
    long double alpha_r, alpha_i;
    long double safmax = 1./safmin;
    long double sigma;
    long double maxab = MAX(ada,adb);
    long double scale = MIN(MAX(safmin,maxab), safmax);


    long double aa_r    = da_r / scale;
    long double aa_i    = da_i / scale;
    long double bb_r    = db_r / scale;
    long double bb_i    = db_i / scale;

    if (ada > adb)
        sigma = copysign(1.,da_r);
    else
        sigma = copysign(1.,db_r);

    r = sigma * scale * sqrt(aa_r * aa_r + aa_i * aa_i + bb_r * bb_r + bb_i * bb_i);


    alpha_r = da_r / ada;
    alpha_i = da_i / ada;

    *(C + 0)  = ada / r;
    *(S + 0)  = (alpha_r * db_r + alpha_i *db_i) / r;
    *(S + 1)  = (alpha_i * db_r - alpha_r *db_i) / r;
    *(DA + 0) = alpha_r * r;
    *(DA + 1) = alpha_i * r;
  }
#else
  FLOAT da_r = *(DA + 0);
  FLOAT da_i = *(DA + 1);
  FLOAT db_r = *(DB + 0);
  FLOAT db_i = *(DB + 1);
  FLOAT r;

  FLOAT ada = fabs(da_r) + fabs(da_i);
  FLOAT ada = fabs(db_r) + fabs(db_i);

  PRINT_DEBUG_NAME;

  IDEBUG_START;

  FUNCTION_PROFILE_START();

  if (ada == ZERO) {
    *C        = ZERO;
    *(S  + 0) = ONE;
    *(S  + 1) = ZERO;
    *(DA + 0) = db_r;
    *(DA + 1) = db_i;
  } else {
    long double safmax = 1./safmin;
    FLOAT scale;
    FLOAT aa_r, aa_i, bb_r, bb_i;
    FLOAT alpha_r, alpha_i;

    aa_r = fabs(da_r);
    aa_i = fabs(da_i);

    if (aa_i > aa_r) {
      aa_r = fabs(da_i);
      aa_i = fabs(da_r);
    }

    if (aa_r == ZERO) {
	ada = 0.;
    } else {
        scale = (aa_i / aa_r);
        ada = aa_r * sqrt(ONE + scale * scale);
    }

    bb_r = fabs(db_r);
    bb_i = fabs(db_i);

    if (bb_i > bb_r) {
      bb_r = fabs(bb_i);
      bb_i = fabs(bb_r);
    }

    if (bb_r == ZERO) {
	adb = 0.;
    } else {
    	scale = (bb_i / bb_r);
    	adb = bb_r * sqrt(ONE + scale * scale);
    }
    FLOAT maxab = MAX(ada,adb);
    scale = MIN(MAX(safmin,maxab), safmax);

    aa_r    = da_r / scale;
    aa_i    = da_i / scale;
    bb_r    = db_r / scale;
    bb_i    = db_i / scale;

    if (ada > adb)
        sigma = copysign(1.,da_r);
    else
        sigma = copysign(1.,db_r);

    r = sigma * scale * sqrt(aa_r * aa_r + aa_i * aa_i + bb_r * bb_r + bb_i * bb_i);

    alpha_r = da_r / ada;
    alpha_i = da_i / ada;

    *(C + 0)  = ada / r;
    *(S + 0)  = (alpha_r * db_r + alpha_i *db_i) / r;
    *(S + 1)  = (alpha_i * db_r - alpha_r *db_i) / r;
    *(DA + 0) = alpha_r * r;
    *(DA + 1) = alpha_i * r;
  }
#endif

  FUNCTION_PROFILE_END(4, 4, 4);

  IDEBUG_END;

  return;
}