/* nag_dgbbrd (f08lec) Example Program. * * Copyright 2001 Numerical Algorithms Group. * * Mark 7, 2001. */ #include #include #include #include #include #include #include #include int main(void) { /* Scalars */ double alpha, beta, norm; Integer i, j, kl, ku, m, n, ncc, pdab, pdaw, pdc, pdf, pdq; Integer pdpt, d_len, e_len; Integer exit_status = 0; NagError fail; Nag_OrderType order; /* Arrays */ double *ab = 0, *aw = 0, *c = 0, *d = 0, *e = 0, *f = 0, *pt = 0, *q = 0; #ifdef NAG_COLUMN_MAJOR #define AB(I, J) ab[(J - 1) * pdab + ku + I - J] #define AW(I, J) aw[(J - 1) * pdaw + I - 1] #define F(I, J) f[(J - 1) * pdf + I - 1] #define Q(I, J) q[(J - 1) * pdq + I - 1] order = Nag_ColMajor; #else #define AB(I, J) ab[(I - 1) * pdab + kl + J - I] #define AW(I, J) aw[(I - 1) * pdaw + J - 1] #define F(I, J) f[(I - 1) * pdf + J - 1] #define Q(I, J) q[(I - 1) * pdq + J - 1] order = Nag_RowMajor; #endif INIT_FAIL(fail); printf("nag_dgbbrd (f08lec) Example Program Results\n"); /* Skip heading in data file */ scanf("%*[^\n] "); scanf("%ld%ld%ld%ld%ld%*[^\n] ", &m, &n, &kl, &ku, &ncc); #ifdef NAG_COLUMN_MAJOR pdab = kl + ku + 1; pdaw = m; pdf = m; pdq = m; pdpt = n; pdc = m; #else pdab = kl + ku + 1; pdaw = n; pdf = n; pdq = m; pdpt = n; pdc = MAX(1, ncc); #endif d_len = MIN(m, n); e_len = MIN(m, n) - 1; /* Allocate memory */ if (!(ab = NAG_ALLOC((kl+ku+1) * m, double)) || !(aw = NAG_ALLOC(m * n, double)) || !(f = NAG_ALLOC(m * n, double)) || !(c = NAG_ALLOC(m * MAX(1, ncc), double)) || !(d = NAG_ALLOC(d_len, double)) || !(e = NAG_ALLOC(e_len, double)) || !(pt = NAG_ALLOC(n * n, double)) || !(q = NAG_ALLOC(m * m, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read A from data file */ for (i = 1; i <= m; ++i) { for (j = MAX(1, i - kl); j <= MIN(n, i + ku); ++j) scanf("%lf", &AB(i, j)); } scanf("%*[^\n] "); /* Copy AB into AW */ for (i = 1; i <= m; ++i) { for (j = 1; j <= n; ++j) { if(j >= MAX(1, i - kl) && j<= MIN(n, i + ku)) AW(i, j) = AB(i, j); else AW(i, j) = 0; } } /* nag_gen_real_mat_print (x04cac): Print Matrix A. */ fflush(stdout); nag_gen_real_mat_print(order, Nag_GeneralMatrix, Nag_NonUnitDiag, m, n, aw, pdaw, "Matrix A", 0, &fail); printf("\n"); if (fail.code != NE_NOERROR) { printf("Error from nag_gen_real_mat_print (x04cac).\n%s\n", fail.message); exit_status = 1; goto END; } /* Reduce A to bidiagonal form */ /* nag_dgbbrd (f08lec). * Reduction of real rectangular band matrix to upper * bidiagonal form */ nag_dgbbrd(order, Nag_FormBoth, m, n, ncc, kl, ku, ab, pdab, d, e, q, pdq, pt, pdpt, c, pdc, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dgbbrd (f08lec).\n%s\n", fail.message); exit_status = 1; goto END; } /* F = Q*B */ for(i = 1; i <= m; i++) { F(i, 1) = Q(i, 1) * d[0]; for(j = 2; j <= n; j++) { F(i, j) = (Q(i, j) * d[j-1]) + (Q(i, j-1) * e[j-2]); } } /* nag_dgemm (f16yac): Compute A - Q*B*P^T from the factorization of A */ /* and store in matrix AW */ alpha = -1.0; beta = 1.0; nag_dgemm(order, Nag_NoTrans, Nag_NoTrans, m, n, n, alpha, f, pdf, pt, pdpt, beta, aw, pdaw, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dgemm (f16yac).\n%s\n", fail.message); exit_status = 1; goto END; } /* nag_dge_norm (f16rac): Find norm of matrix AW and print warning if */ /* it is too large*/ nag_dge_norm(order, Nag_OneNorm, m, n, aw, pdaw, &norm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_dge_norm (f16rac).\n%s\n", fail.message); exit_status = 1; goto END; } if (norm > pow(x02ajc(),0.8)) { printf("%s\n%s\n","Norm of A-(Q*B*P^T) is much greater than 0.", "Schur factorization has failed."); } END: if (ab) NAG_FREE(ab); if (aw) NAG_FREE(aw); if (c) NAG_FREE(c); if (d) NAG_FREE(d); if (e) NAG_FREE(e); if (f) NAG_FREE(f); if (pt) NAG_FREE(pt); if (q) NAG_FREE(q); return exit_status; }