--- /dev/null
+#include "clapack.h"
+
+/* Subroutine */ int dlasr_(char *side, char *pivot, char *direct, integer *m,
+ integer *n, doublereal *c__, doublereal *s, doublereal *a, integer *
+ lda)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, i__1, i__2;
+
+ /* Local variables */
+ integer i__, j, info;
+ doublereal temp;
+ extern logical lsame_(char *, char *);
+ doublereal ctemp, stemp;
+ extern /* Subroutine */ int xerbla_(char *, integer *);
+
+
+/* -- LAPACK auxiliary routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+/* .. Array Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* DLASR applies a sequence of plane rotations to a real matrix A, */
+/* from either the left or the right. */
+
+/* When SIDE = 'L', the transformation takes the form */
+
+/* A := P*A */
+
+/* and when SIDE = 'R', the transformation takes the form */
+
+/* A := A*P**T */
+
+/* where P is an orthogonal matrix consisting of a sequence of z plane */
+/* rotations, with z = M when SIDE = 'L' and z = N when SIDE = 'R', */
+/* and P**T is the transpose of P. */
+
+/* When DIRECT = 'F' (Forward sequence), then */
+
+/* P = P(z-1) * ... * P(2) * P(1) */
+
+/* and when DIRECT = 'B' (Backward sequence), then */
+
+/* P = P(1) * P(2) * ... * P(z-1) */
+
+/* where P(k) is a plane rotation matrix defined by the 2-by-2 rotation */
+
+/* R(k) = ( c(k) s(k) ) */
+/* = ( -s(k) c(k) ). */
+
+/* When PIVOT = 'V' (Variable pivot), the rotation is performed */
+/* for the plane (k,k+1), i.e., P(k) has the form */
+
+/* P(k) = ( 1 ) */
+/* ( ... ) */
+/* ( 1 ) */
+/* ( c(k) s(k) ) */
+/* ( -s(k) c(k) ) */
+/* ( 1 ) */
+/* ( ... ) */
+/* ( 1 ) */
+
+/* where R(k) appears as a rank-2 modification to the identity matrix in */
+/* rows and columns k and k+1. */
+
+/* When PIVOT = 'T' (Top pivot), the rotation is performed for the */
+/* plane (1,k+1), so P(k) has the form */
+
+/* P(k) = ( c(k) s(k) ) */
+/* ( 1 ) */
+/* ( ... ) */
+/* ( 1 ) */
+/* ( -s(k) c(k) ) */
+/* ( 1 ) */
+/* ( ... ) */
+/* ( 1 ) */
+
+/* where R(k) appears in rows and columns 1 and k+1. */
+
+/* Similarly, when PIVOT = 'B' (Bottom pivot), the rotation is */
+/* performed for the plane (k,z), giving P(k) the form */
+
+/* P(k) = ( 1 ) */
+/* ( ... ) */
+/* ( 1 ) */
+/* ( c(k) s(k) ) */
+/* ( 1 ) */
+/* ( ... ) */
+/* ( 1 ) */
+/* ( -s(k) c(k) ) */
+
+/* where R(k) appears in rows and columns k and z. The rotations are */
+/* performed without ever forming P(k) explicitly. */
+
+/* Arguments */
+/* ========= */
+
+/* SIDE (input) CHARACTER*1 */
+/* Specifies whether the plane rotation matrix P is applied to */
+/* A on the left or the right. */
+/* = 'L': Left, compute A := P*A */
+/* = 'R': Right, compute A:= A*P**T */
+
+/* PIVOT (input) CHARACTER*1 */
+/* Specifies the plane for which P(k) is a plane rotation */
+/* matrix. */
+/* = 'V': Variable pivot, the plane (k,k+1) */
+/* = 'T': Top pivot, the plane (1,k+1) */
+/* = 'B': Bottom pivot, the plane (k,z) */
+
+/* DIRECT (input) CHARACTER*1 */
+/* Specifies whether P is a forward or backward sequence of */
+/* plane rotations. */
+/* = 'F': Forward, P = P(z-1)*...*P(2)*P(1) */
+/* = 'B': Backward, P = P(1)*P(2)*...*P(z-1) */
+
+/* M (input) INTEGER */
+/* The number of rows of the matrix A. If m <= 1, an immediate */
+/* return is effected. */
+
+/* N (input) INTEGER */
+/* The number of columns of the matrix A. If n <= 1, an */
+/* immediate return is effected. */
+
+/* C (input) DOUBLE PRECISION array, dimension */
+/* (M-1) if SIDE = 'L' */
+/* (N-1) if SIDE = 'R' */
+/* The cosines c(k) of the plane rotations. */
+
+/* S (input) DOUBLE PRECISION array, dimension */
+/* (M-1) if SIDE = 'L' */
+/* (N-1) if SIDE = 'R' */
+/* The sines s(k) of the plane rotations. The 2-by-2 plane */
+/* rotation part of the matrix P(k), R(k), has the form */
+/* R(k) = ( c(k) s(k) ) */
+/* ( -s(k) c(k) ). */
+
+/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
+/* The M-by-N matrix A. On exit, A is overwritten by P*A if */
+/* SIDE = 'R' or by A*P**T if SIDE = 'L'. */
+
+/* LDA (input) INTEGER */
+/* The leading dimension of the array A. LDA >= max(1,M). */
+
+/* ===================================================================== */
+
+/* .. Parameters .. */
+/* .. */
+/* .. Local Scalars .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. External Subroutines .. */
+/* .. */
+/* .. Intrinsic Functions .. */
+/* .. */
+/* .. Executable Statements .. */
+
+/* Test the input parameters */
+
+ /* Parameter adjustments */
+ --c__;
+ --s;
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+
+ /* Function Body */
+ info = 0;
+ if (! (lsame_(side, "L") || lsame_(side, "R"))) {
+ info = 1;
+ } else if (! (lsame_(pivot, "V") || lsame_(pivot,
+ "T") || lsame_(pivot, "B"))) {
+ info = 2;
+ } else if (! (lsame_(direct, "F") || lsame_(direct,
+ "B"))) {
+ info = 3;
+ } else if (*m < 0) {
+ info = 4;
+ } else if (*n < 0) {
+ info = 5;
+ } else if (*lda < max(1,*m)) {
+ info = 9;
+ }
+ if (info != 0) {
+ xerbla_("DLASR ", &info);
+ return 0;
+ }
+
+/* Quick return if possible */
+
+ if (*m == 0 || *n == 0) {
+ return 0;
+ }
+ if (lsame_(side, "L")) {
+
+/* Form P * A */
+
+ if (lsame_(pivot, "V")) {
+ if (lsame_(direct, "F")) {
+ i__1 = *m - 1;
+ for (j = 1; j <= i__1; ++j) {
+ ctemp = c__[j];
+ stemp = s[j];
+ if (ctemp != 1. || stemp != 0.) {
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp = a[j + 1 + i__ * a_dim1];
+ a[j + 1 + i__ * a_dim1] = ctemp * temp - stemp *
+ a[j + i__ * a_dim1];
+ a[j + i__ * a_dim1] = stemp * temp + ctemp * a[j
+ + i__ * a_dim1];
+/* L10: */
+ }
+ }
+/* L20: */
+ }
+ } else if (lsame_(direct, "B")) {
+ for (j = *m - 1; j >= 1; --j) {
+ ctemp = c__[j];
+ stemp = s[j];
+ if (ctemp != 1. || stemp != 0.) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = a[j + 1 + i__ * a_dim1];
+ a[j + 1 + i__ * a_dim1] = ctemp * temp - stemp *
+ a[j + i__ * a_dim1];
+ a[j + i__ * a_dim1] = stemp * temp + ctemp * a[j
+ + i__ * a_dim1];
+/* L30: */
+ }
+ }
+/* L40: */
+ }
+ }
+ } else if (lsame_(pivot, "T")) {
+ if (lsame_(direct, "F")) {
+ i__1 = *m;
+ for (j = 2; j <= i__1; ++j) {
+ ctemp = c__[j - 1];
+ stemp = s[j - 1];
+ if (ctemp != 1. || stemp != 0.) {
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp = a[j + i__ * a_dim1];
+ a[j + i__ * a_dim1] = ctemp * temp - stemp * a[
+ i__ * a_dim1 + 1];
+ a[i__ * a_dim1 + 1] = stemp * temp + ctemp * a[
+ i__ * a_dim1 + 1];
+/* L50: */
+ }
+ }
+/* L60: */
+ }
+ } else if (lsame_(direct, "B")) {
+ for (j = *m; j >= 2; --j) {
+ ctemp = c__[j - 1];
+ stemp = s[j - 1];
+ if (ctemp != 1. || stemp != 0.) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = a[j + i__ * a_dim1];
+ a[j + i__ * a_dim1] = ctemp * temp - stemp * a[
+ i__ * a_dim1 + 1];
+ a[i__ * a_dim1 + 1] = stemp * temp + ctemp * a[
+ i__ * a_dim1 + 1];
+/* L70: */
+ }
+ }
+/* L80: */
+ }
+ }
+ } else if (lsame_(pivot, "B")) {
+ if (lsame_(direct, "F")) {
+ i__1 = *m - 1;
+ for (j = 1; j <= i__1; ++j) {
+ ctemp = c__[j];
+ stemp = s[j];
+ if (ctemp != 1. || stemp != 0.) {
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp = a[j + i__ * a_dim1];
+ a[j + i__ * a_dim1] = stemp * a[*m + i__ * a_dim1]
+ + ctemp * temp;
+ a[*m + i__ * a_dim1] = ctemp * a[*m + i__ *
+ a_dim1] - stemp * temp;
+/* L90: */
+ }
+ }
+/* L100: */
+ }
+ } else if (lsame_(direct, "B")) {
+ for (j = *m - 1; j >= 1; --j) {
+ ctemp = c__[j];
+ stemp = s[j];
+ if (ctemp != 1. || stemp != 0.) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = a[j + i__ * a_dim1];
+ a[j + i__ * a_dim1] = stemp * a[*m + i__ * a_dim1]
+ + ctemp * temp;
+ a[*m + i__ * a_dim1] = ctemp * a[*m + i__ *
+ a_dim1] - stemp * temp;
+/* L110: */
+ }
+ }
+/* L120: */
+ }
+ }
+ }
+ } else if (lsame_(side, "R")) {
+
+/* Form A * P' */
+
+ if (lsame_(pivot, "V")) {
+ if (lsame_(direct, "F")) {
+ i__1 = *n - 1;
+ for (j = 1; j <= i__1; ++j) {
+ ctemp = c__[j];
+ stemp = s[j];
+ if (ctemp != 1. || stemp != 0.) {
+ i__2 = *m;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp = a[i__ + (j + 1) * a_dim1];
+ a[i__ + (j + 1) * a_dim1] = ctemp * temp - stemp *
+ a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = stemp * temp + ctemp * a[
+ i__ + j * a_dim1];
+/* L130: */
+ }
+ }
+/* L140: */
+ }
+ } else if (lsame_(direct, "B")) {
+ for (j = *n - 1; j >= 1; --j) {
+ ctemp = c__[j];
+ stemp = s[j];
+ if (ctemp != 1. || stemp != 0.) {
+ i__1 = *m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = a[i__ + (j + 1) * a_dim1];
+ a[i__ + (j + 1) * a_dim1] = ctemp * temp - stemp *
+ a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = stemp * temp + ctemp * a[
+ i__ + j * a_dim1];
+/* L150: */
+ }
+ }
+/* L160: */
+ }
+ }
+ } else if (lsame_(pivot, "T")) {
+ if (lsame_(direct, "F")) {
+ i__1 = *n;
+ for (j = 2; j <= i__1; ++j) {
+ ctemp = c__[j - 1];
+ stemp = s[j - 1];
+ if (ctemp != 1. || stemp != 0.) {
+ i__2 = *m;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = ctemp * temp - stemp * a[
+ i__ + a_dim1];
+ a[i__ + a_dim1] = stemp * temp + ctemp * a[i__ +
+ a_dim1];
+/* L170: */
+ }
+ }
+/* L180: */
+ }
+ } else if (lsame_(direct, "B")) {
+ for (j = *n; j >= 2; --j) {
+ ctemp = c__[j - 1];
+ stemp = s[j - 1];
+ if (ctemp != 1. || stemp != 0.) {
+ i__1 = *m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = ctemp * temp - stemp * a[
+ i__ + a_dim1];
+ a[i__ + a_dim1] = stemp * temp + ctemp * a[i__ +
+ a_dim1];
+/* L190: */
+ }
+ }
+/* L200: */
+ }
+ }
+ } else if (lsame_(pivot, "B")) {
+ if (lsame_(direct, "F")) {
+ i__1 = *n - 1;
+ for (j = 1; j <= i__1; ++j) {
+ ctemp = c__[j];
+ stemp = s[j];
+ if (ctemp != 1. || stemp != 0.) {
+ i__2 = *m;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = stemp * a[i__ + *n * a_dim1]
+ + ctemp * temp;
+ a[i__ + *n * a_dim1] = ctemp * a[i__ + *n *
+ a_dim1] - stemp * temp;
+/* L210: */
+ }
+ }
+/* L220: */
+ }
+ } else if (lsame_(direct, "B")) {
+ for (j = *n - 1; j >= 1; --j) {
+ ctemp = c__[j];
+ stemp = s[j];
+ if (ctemp != 1. || stemp != 0.) {
+ i__1 = *m;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ temp = a[i__ + j * a_dim1];
+ a[i__ + j * a_dim1] = stemp * a[i__ + *n * a_dim1]
+ + ctemp * temp;
+ a[i__ + *n * a_dim1] = ctemp * a[i__ + *n *
+ a_dim1] - stemp * temp;
+/* L230: */
+ }
+ }
+/* L240: */
+ }
+ }
+ }
+ }
+
+ return 0;
+
+/* End of DLASR */
+
+} /* dlasr_ */