```01:       SUBROUTINE DLADIV( A, B, C, D, P, Q )
02: *
03: *  -- LAPACK auxiliary routine (version 3.2) --
04: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
05: *     November 2006
06: *
07: *     .. Scalar Arguments ..
08:       DOUBLE PRECISION   A, B, C, D, P, Q
09: *     ..
10: *
11: *  Purpose
12: *  =======
13: *
14: *  DLADIV performs complex division in  real arithmetic
15: *
16: *                        a + i*b
17: *             p + i*q = ---------
18: *                        c + i*d
19: *
20: *  The algorithm is due to Robert L. Smith and can be found
21: *  in D. Knuth, The art of Computer Programming, Vol.2, p.195
22: *
23: *  Arguments
24: *  =========
25: *
26: *  A       (input) DOUBLE PRECISION
27: *  B       (input) DOUBLE PRECISION
28: *  C       (input) DOUBLE PRECISION
29: *  D       (input) DOUBLE PRECISION
30: *          The scalars a, b, c, and d in the above expression.
31: *
32: *  P       (output) DOUBLE PRECISION
33: *  Q       (output) DOUBLE PRECISION
34: *          The scalars p and q in the above expression.
35: *
36: *  =====================================================================
37: *
38: *     .. Local Scalars ..
39:       DOUBLE PRECISION   E, F
40: *     ..
41: *     .. Intrinsic Functions ..
42:       INTRINSIC          ABS
43: *     ..
44: *     .. Executable Statements ..
45: *
46:       IF( ABS( D ).LT.ABS( C ) ) THEN
47:          E = D / C
48:          F = C + D*E
49:          P = ( A+B*E ) / F
50:          Q = ( B-A*E ) / F
51:       ELSE
52:          E = C / D
53:          F = D + C*E
54:          P = ( B+A*E ) / F
55:          Q = ( -A+B*E ) / F
56:       END IF
57: *
58:       RETURN
59: *