1 | ///////////////////////////////////////////////////////////////////////////////// |
---|
2 | //// |
---|
3 | //// Visualization of the nonzero pattern of JtJ in EPS format. |
---|
4 | //// This is part of the sba package, |
---|
5 | //// Copyright (C) 2004 Manolis Lourakis (lourakis@ics.forth.gr) |
---|
6 | //// Institute of Computer Science, Foundation for Research & Technology - Hellas |
---|
7 | //// Heraklion, Crete, Greece. |
---|
8 | //// |
---|
9 | //// This program is free software; you can redistribute it and/or modify |
---|
10 | //// it under the terms of the GNU General Public License as published by |
---|
11 | //// the Free Software Foundation; either version 2 of the License, or |
---|
12 | //// (at your option) any later version. |
---|
13 | //// |
---|
14 | //// This program is distributed in the hope that it will be useful, |
---|
15 | //// but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
16 | //// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
17 | //// GNU General Public License for more details. |
---|
18 | //// |
---|
19 | /////////////////////////////////////////////////////////////////////////////////// |
---|
20 | |
---|
21 | #include <stdio.h> |
---|
22 | #include <stdlib.h> |
---|
23 | #include <math.h> |
---|
24 | #include <time.h> |
---|
25 | |
---|
26 | #include "sba.h" |
---|
27 | |
---|
28 | #if 0 |
---|
29 | #define HORZ_OFFSET 3 // horizontal offset from left edge of paper in PostScript points (~ 1 mm) |
---|
30 | #define VERT_OFFSET 6 // vertical offset from lower edge of paper in PostScript points (~ 2 mm) |
---|
31 | #else |
---|
32 | #define HORZ_OFFSET 0 |
---|
33 | #define VERT_OFFSET 0 |
---|
34 | #endif |
---|
35 | |
---|
36 | #define A4_WIDTH 595 // A4 paper width in PostScript points (210 mm, 1 pt = 25.4/72 mm) |
---|
37 | //#define A4_HEIGHT 842 // A4 paper height in PostScript points (297 mm) |
---|
38 | |
---|
39 | /* 1.1x1.1 looks better than 1.0x1.0 */ |
---|
40 | #define DOT_WIDTH 1.1 |
---|
41 | #define DOT_HEIGHT 1.1 |
---|
42 | |
---|
43 | /* saves the nonzero pattern of JtJ in EPS format. If mcon!=0, then all U_j and W_ij with j<mcon |
---|
44 | * are assumed to be zero and are therefore clipped off. |
---|
45 | */ |
---|
46 | void sba_hessian2eps(int n, int m, int mcon, int cnp, int pnp, struct sba_crsm *idxij, int *rcidxs, int *rcsubs, |
---|
47 | double *U, double *V, double *W, char *fname) |
---|
48 | { |
---|
49 | register int i, j, ii, jj; |
---|
50 | int i0, j0; |
---|
51 | int mmcon, nvars, Usz, Vsz, Wsz, nnz, bbox=1, denseblocks=0; |
---|
52 | register double *ptr1, *ptr2; |
---|
53 | FILE *fp; |
---|
54 | time_t tim; |
---|
55 | |
---|
56 | mmcon=m-mcon; |
---|
57 | nvars=mmcon*cnp + n*pnp; |
---|
58 | Usz=cnp * cnp; Vsz=pnp * pnp; Wsz=cnp * pnp; |
---|
59 | |
---|
60 | if((fp=fopen(fname, "w"))==NULL){ |
---|
61 | fprintf(stderr, "JtJ2EPS(): failed to open file %s for writing\n", fname); |
---|
62 | exit(1); |
---|
63 | } |
---|
64 | |
---|
65 | /* print EPS preamble */ |
---|
66 | time(&tim); |
---|
67 | fprintf(fp, "%%!PS-Adobe-2.0 EPSF-2.0\n%%%%Title: %s\n%%%%Creator: sba ver. %s\n", fname, SBA_VERSION); |
---|
68 | fprintf(fp, "%%%%CreationDate: %s", ctime(&tim)); |
---|
69 | fprintf(fp, "%%%%BoundingBox: 0 0 %d %d\n%%%%Magnification: 1.0000\n%%%%Page: 1 1\n%%%%EndComments\n\n", A4_WIDTH, A4_WIDTH+VERT_OFFSET); |
---|
70 | fprintf(fp, "/origstate save def\ngsave\n0 setgray\n"); |
---|
71 | if(bbox) fprintf(fp, "0 0 %d %d rectstroke\n", A4_WIDTH, A4_WIDTH+VERT_OFFSET); |
---|
72 | |
---|
73 | /* move the coordinate system to the upper left corner with axes pointing as shown below: |
---|
74 | * +------> y |
---|
75 | * | |
---|
76 | * | |
---|
77 | * | |
---|
78 | * v x |
---|
79 | */ |
---|
80 | fprintf(fp, "%d %d translate\n-90 rotate\n", HORZ_OFFSET, A4_WIDTH+VERT_OFFSET); |
---|
81 | /* scale the coordinate system so that the hessian pattern fits in the narrowest page dimension */ |
---|
82 | fprintf(fp, "%.4lf %.4lf scale\n", ((double)(A4_WIDTH-HORZ_OFFSET))/nvars, ((double)(A4_WIDTH-HORZ_OFFSET))/nvars); |
---|
83 | |
---|
84 | /* define dot dimensions */ |
---|
85 | fprintf(fp, "/w %g def\n/h %g def\n", DOT_WIDTH, DOT_HEIGHT); |
---|
86 | /* define shorthand for rectfill */ |
---|
87 | fprintf(fp, "/R { rectfill } bind def\n\n"); |
---|
88 | |
---|
89 | /* process block row j: [0, ..., 0, U_j, 0, ..., 0, W1j, ..., Wnj] */ |
---|
90 | for(j=mcon; j<m; ++j){ |
---|
91 | ptr1=U + j*Usz; // set ptr1 to point to U_j |
---|
92 | nnz=sba_crsm_col_elmidxs(idxij, j, rcidxs, rcsubs); /* find nonzero W_ij, i=0...n-1 */ |
---|
93 | |
---|
94 | i0=(j-mcon)*cnp; |
---|
95 | for(ii=0; ii<cnp; ++ii){ |
---|
96 | j0=(j-mcon)*cnp; |
---|
97 | for(jj=0; jj<cnp; ++jj){ // row ii of U_j |
---|
98 | if(denseblocks || (!denseblocks && ptr1[ii*cnp+jj]!=0.0)) |
---|
99 | fprintf(fp, "%d %d w h R\n", i0+ii, j0+jj); |
---|
100 | } |
---|
101 | |
---|
102 | for(i=0; i<nnz; ++i){ |
---|
103 | /* set ptr2 to point to W_ij, actual row number in rcsubs[i] */ |
---|
104 | ptr2=W + idxij->val[rcidxs[i]]*Wsz; |
---|
105 | j0=mmcon*cnp+rcsubs[i]*pnp; |
---|
106 | for(jj=0; jj<pnp; ++jj){ // row ii of W_ij |
---|
107 | if(denseblocks || (!denseblocks && ptr2[ii*pnp+jj]!=0.0)) |
---|
108 | fprintf(fp, "%d %d w h R\n", i0+ii, j0+jj); |
---|
109 | } |
---|
110 | } |
---|
111 | } |
---|
112 | } |
---|
113 | |
---|
114 | /* process block row mmcon+i: [W_i1^T, ..., W_im^T, 0, ..., 0, V_i, 0, ..., 0] */ |
---|
115 | for(i=0; i<n; ++i){ |
---|
116 | nnz=sba_crsm_row_elmidxs(idxij, i, rcidxs, rcsubs); /* find nonzero W_ij, j=0...m-1 */ |
---|
117 | ptr1=V + i*Vsz; // set ptr1 to point to V_i |
---|
118 | |
---|
119 | i0=cnp*mmcon+i*pnp; |
---|
120 | for(ii=0; ii<pnp; ++ii){ |
---|
121 | for(j=0; j<nnz; ++j){ |
---|
122 | /* set ptr2 to point to W_ij, actual column number in rcsubs[j] */ |
---|
123 | if(rcsubs[j]<mcon) continue; /* W_ij is zero */ |
---|
124 | |
---|
125 | ptr2=W + idxij->val[rcidxs[j]]*Wsz; |
---|
126 | |
---|
127 | j0=(rcsubs[j]-mcon)*cnp; |
---|
128 | for(jj=0; jj<cnp; ++jj){ // row ii of W_ij^T |
---|
129 | if(denseblocks || (!denseblocks && ptr2[jj*pnp+ii]!=0.0)) |
---|
130 | fprintf(fp, "%d %d w h R\n", i0+ii, j0+jj); |
---|
131 | } |
---|
132 | } |
---|
133 | |
---|
134 | j0=cnp*mmcon+i*pnp; |
---|
135 | for(jj=0; jj<pnp; ++jj){ // row ii of V_i |
---|
136 | if(denseblocks || (!denseblocks && ptr1[ii*pnp+jj]!=0.0)) |
---|
137 | fprintf(fp, "%d %d w h R\n", i0+ii, j0+jj); |
---|
138 | } |
---|
139 | } |
---|
140 | } |
---|
141 | |
---|
142 | fprintf(fp, "grestore\norigstate restore\n\n%%%%Trailer"); |
---|
143 | fclose(fp); |
---|
144 | } |
---|