Mercurial > octave-libgccjit
view libcruft/amos/zwrsk.f @ 8964:f4f4d65faaa0
Implement sparse * diagonal and diagonal * sparse operations, double-prec only.
Date: Sun, 8 Mar 2009 16:28:18 -0400
These preserve sparsity, so eye(5) * sprand (5, 5, .2) is *sparse*
and not dense. This may affect people who use multiplication by
eye() rather than full().
The liboctave routines do *not* check if arguments are scalars in
disguise. There is a type problem with checking at that level. I
suspect we want diag * "sparse scalar" to stay diagonal, but we have
to return a sparse matrix at the liboctave. Rather than worrying
about that in liboctave, we cope with it when binding to Octave and
return the correct higher-level type.
The implementation is in Sparse-diag-op-defs.h rather than
Sparse-op-defs.h to limit recompilation. And the implementations
are templates rather than macros to produce better compiler errors
and debugging information.
| author | Jason Riedy <jason@acm.org> |
|---|---|
| date | Mon, 09 Mar 2009 17:49:13 -0400 |
| parents | 8b0cb8f79fdc |
| children |
line wrap: on
line source
SUBROUTINE ZWRSK(ZRR, ZRI, FNU, KODE, N, YR, YI, NZ, CWR, CWI, * TOL, ELIM, ALIM) C***BEGIN PROLOGUE ZWRSK C***REFER TO ZBESI,ZBESK C C ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY C NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN C C***ROUTINES CALLED D1MACH,ZBKNU,ZRATI,XZABS C***END PROLOGUE ZWRSK C COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR DOUBLE PRECISION ACT, ACW, ALIM, ASCLE, CINUI, CINUR, CSCLR, CTI, * CTR, CWI, CWR, C1I, C1R, C2I, C2R, ELIM, FNU, PTI, PTR, RACT, * STI, STR, TOL, YI, YR, ZRI, ZRR, XZABS, D1MACH INTEGER I, KODE, N, NW, NZ DIMENSION YR(N), YI(N), CWR(2), CWI(2) C----------------------------------------------------------------------- C I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS C Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE C WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU. C----------------------------------------------------------------------- NZ = 0 CALL ZBKNU(ZRR, ZRI, FNU, KODE, 2, CWR, CWI, NW, TOL, ELIM, ALIM) IF (NW.NE.0) GO TO 50 CALL ZRATI(ZRR, ZRI, FNU, N, YR, YI, TOL) C----------------------------------------------------------------------- C RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z), C R(FNU+J-1,Z)=Y(J), J=1,...,N C----------------------------------------------------------------------- CINUR = 1.0D0 CINUI = 0.0D0 IF (KODE.EQ.1) GO TO 10 CINUR = DCOS(ZRI) CINUI = DSIN(ZRI) 10 CONTINUE C----------------------------------------------------------------------- C ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH C THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE C SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT C THE RESULT IS ON SCALE. C----------------------------------------------------------------------- ACW = XZABS(CWR(2),CWI(2)) ASCLE = 1.0D+3*D1MACH(1)/TOL CSCLR = 1.0D0 IF (ACW.GT.ASCLE) GO TO 20 CSCLR = 1.0D0/TOL GO TO 30 20 CONTINUE ASCLE = 1.0D0/ASCLE IF (ACW.LT.ASCLE) GO TO 30 CSCLR = TOL 30 CONTINUE C1R = CWR(1)*CSCLR C1I = CWI(1)*CSCLR C2R = CWR(2)*CSCLR C2I = CWI(2)*CSCLR STR = YR(1) STI = YI(1) C----------------------------------------------------------------------- C CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0D0/CABS(CT) PREVENTS C UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT) C----------------------------------------------------------------------- PTR = STR*C1R - STI*C1I PTI = STR*C1I + STI*C1R PTR = PTR + C2R PTI = PTI + C2I CTR = ZRR*PTR - ZRI*PTI CTI = ZRR*PTI + ZRI*PTR ACT = XZABS(CTR,CTI) RACT = 1.0D0/ACT CTR = CTR*RACT CTI = -CTI*RACT PTR = CINUR*RACT PTI = CINUI*RACT CINUR = PTR*CTR - PTI*CTI CINUI = PTR*CTI + PTI*CTR YR(1) = CINUR*CSCLR YI(1) = CINUI*CSCLR IF (N.EQ.1) RETURN DO 40 I=2,N PTR = STR*CINUR - STI*CINUI CINUI = STR*CINUI + STI*CINUR CINUR = PTR STR = YR(I) STI = YI(I) YR(I) = CINUR*CSCLR YI(I) = CINUI*CSCLR 40 CONTINUE RETURN 50 CONTINUE NZ = -1 IF(NW.EQ.(-2)) NZ=-2 RETURN END