/* nag_sparse_nsym_fac_solve(f11dcc) Example Program. * * Copyright 1998 Numerical Algorithms Group. * * Mark 5, 1998. * */ #include #include #include #include #include int main(void) { double dtol; double *a=0, *b=0; double *x=0; double rnorm; double tol; Integer *irow, *icol; Integer *istr=0, *idiag, *ipivp=0, *ipivq=0; Integer i, m, n, nnzc; Integer lfill, npivm; Integer maxitn; Integer itn; Integer nnz; Integer num; Nag_SparseNsym_Method method; Nag_SparseNsym_Piv pstrat; Nag_SparseNsym_Fact milu; Nag_Sparse_Comm comm; char char_enum[20]; Vprintf("f11dcc Example Program Results\n"); /* Skip heading in data file */ Vscanf("%*[^\n]"); Vscanf("%ld%*[^\n]",&n); Vscanf("%ld%*[^\n]",&nnz); Vscanf("%s%*[^\n]",char_enum); if (!strcmp(char_enum, "RGMRES")) method = Nag_SparseNsym_RGMRES; else if (!strcmp(char_enum, "CGS")) method = Nag_SparseNsym_CGS; else if (!strcmp(char_enum, "BiCGSTAB")) method = Nag_SparseNsym_BiCGSTAB; else { Vprintf("Unrecognised string for method enum representation.\n"); return EXIT_FAILURE; } Vscanf("%ld%lf%*[^\n]",&lfill,&dtol); Vscanf("%s%*[^\n]", char_enum); if (!strcmp(char_enum, "NoPiv")) pstrat = Nag_SparseNsym_NoPiv; else if (!strcmp(char_enum, "UserPiv")) pstrat = Nag_SparseNsym_UserPiv; else if (!strcmp(char_enum, "PartialPiv")) pstrat = Nag_SparseNsym_PartialPiv; else if (!strcmp(char_enum, "CompletePiv")) pstrat = Nag_SparseNsym_CompletePiv; else { Vprintf("Unrecognised string for psstrat enum representation.\n"); return EXIT_FAILURE; } Vscanf("%s%*[^\n]", char_enum); if (!strcmp(char_enum, "ModFact")) milu = Nag_SparseNsym_ModFact; else if (!strcmp(char_enum, "UnModFact")) milu = Nag_SparseNsym_UnModFact; else { Vprintf("Unrecognised string for milu enum representation.\n"); return EXIT_FAILURE; } Vscanf("%ld%lf%ld%*[^\n]",&m,&tol,&maxitn); /* Read the matrix a */ num = 2*nnz; istr = NAG_ALLOC(n+1, Integer); idiag = NAG_ALLOC(n, Integer); ipivp = NAG_ALLOC(n, Integer); ipivq = NAG_ALLOC(n, Integer); x = NAG_ALLOC(n,double); b = NAG_ALLOC(n,double); a = NAG_ALLOC(num,double); irow = NAG_ALLOC(num,Integer); icol = NAG_ALLOC(num,Integer); if (!istr || !idiag || !ipivp || !ipivq ||!irow || !icol || !a || !x || !b) { Vprintf("Allocation failure\n"); return EXIT_FAILURE; } for (i = 1; i <= nnz; ++i) Vscanf("%lf%ld%ld%*[^\n]",&a[i-1], &irow[i-1], &icol[i-1]); /* Read right-hand side vector b and initial approximate solution x */ for (i = 1; i <= n; ++i) Vscanf("%lf",&b[i-1]); Vscanf("%*[^\n]"); for (i = 1; i <= n; ++i) Vscanf("%lf",&x[i-1]); Vscanf("%*[^\n]"); /* Calculate incomplete LU factorization */ f11dac(n, nnz, &a, &num, &irow, &icol, lfill, dtol, pstrat, milu, ipivp, ipivq, istr, idiag, &nnzc, &npivm, NAGERR_DEFAULT); /* Solve Ax = b using F11DCC */ f11dcc(method, n, nnz, a, num, irow, icol, ipivp, ipivq, istr, idiag, b, m, tol, maxitn, x, &rnorm, &itn, &comm, NAGERR_DEFAULT); Vprintf("%s%10ld%s\n","Converged in",itn," iterations"); Vprintf("%s%16.3e\n","Final residual norm =",rnorm); /* Output x */ Vprintf(" x\n"); for (i = 1; i <= n; ++i) Vprintf(" %16.6e\n",x[i-1]); NAG_FREE(istr); NAG_FREE(idiag); NAG_FREE(ipivp); NAG_FREE(ipivq); NAG_FREE(irow); NAG_FREE(icol); NAG_FREE(a); NAG_FREE(x); NAG_FREE(b); return EXIT_SUCCESS; }