NEURON
spoutput.cpp
Go to the documentation of this file.
1 /*
2  * MATRIX OUTPUT MODULE
3  *
4  * Author: Advisor:
5  * Kenneth S. Kundert Alberto Sangiovanni-Vincentelli
6  * UC Berkeley
7  *
8  * This file contains the output-to-file and output-to-screen routines for
9  * the matrix package.
10  *
11  * >>> User accessible functions contained in this file:
12  * spPrint
13  * spFileMatrix
14  * spFileVector
15  * spFileStats
16  *
17  * >>> Other functions contained in this file:
18  */
19 
20 /*
21  * Revision and copyright information.
22  *
23  * Copyright (c) 1985,86,87,88
24  * by Kenneth S. Kundert and the University of California.
25  *
26  * Permission to use, copy, modify, and distribute this software and
27  * its documentation for any purpose and without fee is hereby granted,
28  * provided that the copyright notices appear in all copies and
29  * supporting documentation and that the authors and the University of
30  * California are properly credited. The authors and the University of
31  * California make no representations as to the suitability of this
32  * software for any purpose. It is provided `as is', without express
33  * or implied warranty.
34  */
35 
36 /*
37  * IMPORTS
38  *
39  * >>> Import descriptions:
40  * spconfig.h
41  * Macros that customize the sparse matrix routines.
42  * spmatrix.h
43  * Macros and declarations to be imported by the user.
44  * spdefs.h
45  * Matrix type and macro definitions for the sparse matrix routines.
46  */
47 
48 #define spINSIDE_SPARSE
49 #include "spconfig.h"
50 #include "spdefs.h"
51 #include "spmatrix.h"
52 #include <cfloat>
53 
54 #if DOCUMENTATION
55 
56 /*
57  * PRINT MATRIX
58  *
59  * Formats and send the matrix to standard output. Some elementary
60  * statistics are also output. The matrix is output in a format that is
61  * readable by people.
62  *
63  * >>> Arguments:
64  * Matrix <input> (char *)
65  * Pointer to matrix.
66  * PrintReordered <input> (int)
67  * Indicates whether the matrix should be printed out in its original
68  * form, as input by the user, or whether it should be printed in its
69  * reordered form, as used by the matrix routines. A zero indicates that
70  * the matrix should be printed as inputed, a one indicates that it
71  * should be printed reordered.
72  * Data <input> (int)
73  * Boolean flag that when false indicates that output should be
74  * compressed such that only the existence of an element should be
75  * indicated rather than giving the actual value. Thus 11 times as
76  * many can be printed on a row. A zero signifies that the matrix
77  * should be printed compressed. A one indicates that the matrix
78  * should be printed in all its glory.
79  * Header <input> (int)
80  * Flag indicating that extra information should be given, such as row
81  * and column numbers.
82  *
83  * >>> Local variables:
84  * Col (int)
85  * Column being printed.
86  * ElementCount (int)
87  * Variable used to count the number of nonzero elements in the matrix.
88  * LargestElement (RealNumber)
89  * The magnitude of the largest element in the matrix.
90  * LargestDiag (RealNumber)
91  * The magnitude of the largest diagonal in the matrix.
92  * Magnitude (RealNumber)
93  * The absolute value of the matrix element being printed.
94  * PrintOrdToIntColMap (int [])
95  * A translation array that maps the order that columns will be
96  * printed in (if not PrintReordered) to the internal column numbers.
97  * PrintOrdToIntRowMap (int [])
98  * A translation array that maps the order that rows will be
99  * printed in (if not PrintReordered) to the internal row numbers.
100  * pElement (ElementPtr)
101  * Pointer to the element in the matrix that is to be printed.
102  * pImagElements (ElementPtr [ ])
103  * Array of pointers to elements in the matrix. These pointers point
104  * to the elements whose real values have just been printed. They are
105  * used to quickly access those same elements so their imaginary values
106  * can be printed.
107  * Row (int)
108  * Row being printed.
109  * Size (int)
110  * The size of the matrix.
111  * SmallestDiag (RealNumber)
112  * The magnitude of the smallest diagonal in the matrix.
113  * SmallestElement (RealNumber)
114  * The magnitude of the smallest element in the matrix excluding zero
115  * elements.
116  * StartCol (int)
117  * The column number of the first column to be printed in the group of
118  * columns currently being printed.
119  * StopCol (int)
120  * The column number of the last column to be printed in the group of
121  * columns currently being printed.
122  * Top (int)
123  * The largest expected external row or column number.
124  */
125 
126 void spPrint(char* eMatrix, int PrintReordered, int Data, int Header)
127 {
128  MatrixPtr Matrix = (MatrixPtr)eMatrix;
129  int J = 0;
130  int I, Row, Col, Size, Top, StartCol = 1, StopCol, Columns, ElementCount = 0;
131  double Magnitude, SmallestDiag, SmallestElement;
132  double LargestElement = 0.0, LargestDiag = 0.0;
133  ElementPtr pElement, pImagElements[PRINTER_WIDTH / 10 + 1];
134  int *PrintOrdToIntRowMap, *PrintOrdToIntColMap;
135 
136  /* Begin `spPrint'. */
138  Size = Matrix->Size;
139 
140 /* Create a packed external to internal row and column translation array. */
141  Top = Matrix->AllocatedSize;
142  CALLOC(PrintOrdToIntRowMap, int, Top + 1);
143  CALLOC(PrintOrdToIntColMap, int, Top + 1);
144  if (PrintOrdToIntRowMap == NULL OR PrintOrdToIntColMap == NULL) {
145  Matrix->Error = spNO_MEMORY;
146  return;
147  }
148  for (I = 1; I <= Size; I++) {
149  PrintOrdToIntRowMap[Matrix->IntToExtRowMap[I]] = I;
150  PrintOrdToIntColMap[Matrix->IntToExtColMap[I]] = I;
151  }
152 
153  /* Pack the arrays. */
154  for (J = 1, I = 1; I <= Top; I++) {
155  if (PrintOrdToIntRowMap[I] != 0)
156  PrintOrdToIntRowMap[J++] = PrintOrdToIntRowMap[I];
157  }
158  for (J = 1, I = 1; I <= Top; I++) {
159  if (PrintOrdToIntColMap[I] != 0)
160  PrintOrdToIntColMap[J++] = PrintOrdToIntColMap[I];
161  }
162 
163  /* Print header. */
164  if (Header) {
165  printf("MATRIX SUMMARY\n\n");
166  printf("Size of matrix = %1u x %1u.\n", Size, Size);
167  if (Matrix->Reordered AND PrintReordered)
168  printf("Matrix has been reordered.\n");
169  printf("\n");
170 
171  if (Matrix->Factored)
172  printf("Matrix after factorization:\n");
173  else
174  printf("Matrix before factorization:\n");
175 
176  SmallestElement = DBL_MAX;
177  SmallestDiag = SmallestElement;
178  }
179 
180  /* Determine how many columns to use. */
181  Columns = PRINTER_WIDTH;
182  if (Header)
183  Columns -= 5;
184  if (Data)
185  Columns = (Columns + 1) / 10;
186 
187  /*
188  * Print matrix by printing groups of complete columns until all the columns
189  * are printed.
190  */
191  J = 0;
192  while (J <= Size)
193 
194  /* Calculate index of last column to printed in this group. */
195  {
196  StopCol = StartCol + Columns - 1;
197  if (StopCol > Size)
198  StopCol = Size;
199 
200  /* Label the columns. */
201  if (Header) {
202  if (Data) {
203  printf(" ");
204  for (I = StartCol; I <= StopCol; I++) {
205  if (PrintReordered)
206  Col = I;
207  else
208  Col = PrintOrdToIntColMap[I];
209  printf(" %9d", Matrix->IntToExtColMap[Col]);
210  }
211  printf("\n\n");
212  } else {
213  if (PrintReordered)
214  printf("Columns %1d to %1d.\n", StartCol, StopCol);
215  else {
216  printf("Columns %1d to %1d.\n",
217  Matrix->IntToExtColMap[PrintOrdToIntColMap[StartCol]],
218  Matrix->IntToExtColMap[PrintOrdToIntColMap[StopCol]]);
219  }
220  }
221  }
222 
223  /* Print every row ... */
224  for (I = 1; I <= Size; I++) {
225  if (PrintReordered)
226  Row = I;
227  else
228  Row = PrintOrdToIntRowMap[I];
229 
230  if (Header) {
231  if (PrintReordered AND NOT Data)
232  printf("%4d", I);
233  else
234  printf("%4d", Matrix->IntToExtRowMap[Row]);
235  if (NOT Data)
236  printf(" ");
237  }
238 
239  /* ... in each column of the group. */
240  for (J = StartCol; J <= StopCol; J++) {
241  if (PrintReordered)
242  Col = J;
243  else
244  Col = PrintOrdToIntColMap[J];
245 
246  pElement = Matrix->FirstInCol[Col];
247  while (pElement != NULL AND pElement->Row != Row)
248  pElement = pElement->NextInCol;
249 
250  if (Data)
251  pImagElements[J - StartCol] = pElement;
252 
253  if (pElement != NULL)
254 
255  /* Case where element exists */
256  {
257  if (Data)
258  printf(" %9.3lg", (double)pElement->Real);
259  else
260  printf("x");
261 
262  /* Update status variables */
263  if ((Magnitude = ELEMENT_MAG(pElement)) > LargestElement)
264  LargestElement = Magnitude;
265  if ((Magnitude < SmallestElement) AND(Magnitude != 0.0))
266  SmallestElement = Magnitude;
267  ElementCount++;
268  }
269 
270  /* Case where element is structurally zero */
271  else {
272  if (Data)
273  printf(" ...");
274  else
275  printf(".");
276  }
277  }
278  printf("\n");
279 
280  }
281 
282  /* Calculate index of first column in next group. */
283  StartCol = StopCol;
284  StartCol++;
285  printf("\n");
286  }
287  if (Header) {
288  printf("\nLargest element in matrix = %-1.4lg.\n", LargestElement);
289  printf("Smallest element in matrix = %-1.4lg.\n", SmallestElement);
290 
291  /* Search for largest and smallest diagonal values */
292  for (I = 1; I <= Size; I++) {
293  if (Matrix->Diag[I] != NULL) {
294  Magnitude = ELEMENT_MAG(Matrix->Diag[I]);
295  if (Magnitude > LargestDiag)
296  LargestDiag = Magnitude;
297  if (Magnitude < SmallestDiag)
298  SmallestDiag = Magnitude;
299  }
300  }
301 
302  /* Print the largest and smallest diagonal values */
303  if (Matrix->Factored) {
304  printf("\nLargest diagonal element = %-1.4lg.\n", LargestDiag);
305  printf("Smallest diagonal element = %-1.4lg.\n", SmallestDiag);
306  } else {
307  printf("\nLargest pivot element = %-1.4lg.\n", LargestDiag);
308  printf("Smallest pivot element = %-1.4lg.\n", SmallestDiag);
309  }
310 
311  /* Calculate and print sparsity and number of fill-ins created. */
312  printf("\nDensity = %2.2lf%%.\n", ((double)(ElementCount * 100)) / ((double)(Size * Size)));
313  if (NOT Matrix->NeedsOrdering)
314  printf("Number of fill-ins = %1d.\n", Matrix->Fillins);
315  }
316  printf("\n");
317  (void)fflush(stdout);
318 
319  FREE(PrintOrdToIntColMap);
320  FREE(PrintOrdToIntRowMap);
321  return;
322 }
323 
324 /*
325  * OUTPUT MATRIX TO FILE
326  *
327  * Writes matrix to file in format suitable to be read back in by the
328  * matrix test program.
329  *
330  * >>> Returns:
331  * One is returned if routine was successful, otherwise zero is returned.
332  * The calling function can query errno (the system global error variable)
333  * as to the reason why this routine failed.
334  *
335  * >>> Arguments:
336  * Matrix <input> (char *)
337  * Pointer to matrix.
338  * File <input> (char *)
339  * Name of file into which matrix is to be written.
340  * Label <input> (char *)
341  * String that is transferred to file and is used as a label.
342  * Reordered <input> (BOOLEAN)
343  * Specifies whether matrix should be output in reordered form,
344  * or in original order.
345  * Data <input> (BOOLEAN)
346  * Indicates that the element values should be output along with
347  * the indices for each element. This parameter must be true if
348  * matrix is to be read by the sparse test program.
349  * Header <input> (BOOLEAN)
350  * Indicates that header is desired. This parameter must be true if
351  * matrix is to be read by the sparse test program.
352  *
353  * >>> Local variables:
354  * Col (int)
355  * The original column number of the element being output.
356  * pElement (ElementPtr)
357  * Pointer to an element in the matrix.
358  * pMatrixFile (FILE *)
359  * File pointer to the matrix file.
360  * Row (int)
361  * The original row number of the element being output.
362  * Size (int)
363  * The size of the matrix.
364  */
365 
366 int spFileMatrix(char* eMatrix, char* File, char* Label, int Reordered, int Data, int Header)
367 {
368  MatrixPtr Matrix = (MatrixPtr)eMatrix;
369  int I, Size;
370  ElementPtr pElement;
371  int Row, Col, Err = 0;
372  FILE* pMatrixFile;
373 
374  /* Begin `spFileMatrix'. */
376 
377  /* Open file matrix file in write mode. */
378  if ((pMatrixFile = fopen(File, "w")) == NULL)
379  return 0;
380 
381  /* Output header. */
382  Size = Matrix->Size;
383  if (Header) {
384  if (Matrix->Factored AND Data) {
385  Err = fprintf(pMatrixFile,
386  "Warning : The following matrix is factored in to LU form.\n");
387  }
388  if (Err < 0)
389  return 0;
390  if (fprintf(pMatrixFile, "%s\n", Label) < 0)
391  return 0;
392  Err = fprintf(pMatrixFile, "%d\t%s\n", Size,
393  (Matrix->Complex ? "complex" : "real"));
394  if (Err < 0)
395  return 0;
396  }
397 
398  /* Output matrix. */
399  if (NOT Data) {
400  for (I = 1; I <= Size; I++) {
401  pElement = Matrix->FirstInCol[I];
402  while (pElement != NULL) {
403  if (Reordered) {
404  Row = pElement->Row;
405  Col = I;
406  } else {
407  Row = Matrix->IntToExtRowMap[pElement->Row];
408  Col = Matrix->IntToExtColMap[I];
409  }
410  pElement = pElement->NextInCol;
411  if (fprintf(pMatrixFile, "%d\t%d\n", Row, Col) < 0)
412  return 0;
413  }
414  }
415  /* Output terminator, a line of zeros. */
416  if (Header)
417  if (fprintf(pMatrixFile, "0\t0\n") < 0)
418  return 0;
419  }
420 
421 
422 #if REAL
423  if (Data AND NOT Matrix->Complex) {
424  for (I = 1; I <= Size; I++) {
425  pElement = Matrix->FirstInCol[I];
426  while (pElement != NULL) {
427  Row = Matrix->IntToExtRowMap[pElement->Row];
428  Col = Matrix->IntToExtColMap[I];
429  Err = fprintf(pMatrixFile, "%d\t%d\t%-.15lg\n",
430  Row, Col, (double)pElement->Real);
431  if (Err < 0)
432  return 0;
433  pElement = pElement->NextInCol;
434  }
435  }
436  /* Output terminator, a line of zeros. */
437  if (Header)
438  if (fprintf(pMatrixFile, "0\t0\t0.0\n") < 0)
439  return 0;
440  }
441 #endif /* REAL */
442 
443  /* Close file. */
444  if (fclose(pMatrixFile) < 0)
445  return 0;
446  return 1;
447 }
448 
449 /*
450  * OUTPUT SOURCE VECTOR TO FILE
451  *
452  * Writes vector to file in format suitable to be read back in by the
453  * matrix test program. This routine should be executed after the function
454  * spFileMatrix.
455  *
456  * >>> Returns:
457  * One is returned if routine was successful, otherwise zero is returned.
458  * The calling function can query errno (the system global error variable)
459  * as to the reason why this routine failed.
460  *
461  * >>> Arguments:
462  * Matrix <input> (char *)
463  * Pointer to matrix.
464  * File <input> (char *)
465  * Name of file into which matrix is to be written.
466  * RHS <input> (RealNumber [])
467  * Right-hand side vector.
468  * iRHS <input> (RealNumber [])
469  * Right-hand side vector, imaginary portion.
470  *
471  * >>> Local variables:
472  * pMatrixFile (FILE *)
473  * File pointer to the matrix file.
474  * Size (int)
475  * The size of the matrix.
476  *
477  */
478 
479 int spFileVector(char* eMatrix, char* File, RealVector RHS, RealVector iRHS)
480 {
481  MatrixPtr Matrix = (MatrixPtr)eMatrix;
482  int I, Size;
483  FILE* pMatrixFile;
484 
485  /* Begin `spFileVector'. */
487 
488  /* Open File in append mode. */
489  if ((pMatrixFile = fopen(File, "a")) == NULL)
490  return 0;
491 
492 /* Correct array pointers for ARRAY_OFFSET. */
493 #if NOT ARRAY_OFFSET
494  --RHS;
495 #endif /* NOT ARRAY_OFFSET */
496 
497  /* Output vector. */
498  Size = Matrix->Size;
499 #if REAL
500  {
501  for (I = 1; I <= Size; I++) {
502  if (fprintf(pMatrixFile, "%-.15lg\n", (double)RHS[I]) < 0)
503  return 0;
504  }
505  }
506 #endif /* REAL */
507 
508  /* Close file. */
509  if (fclose(pMatrixFile) < 0)
510  return 0;
511  return 1;
512 }
513 
514 /*
515  * OUTPUT STATISTICS TO FILE
516  *
517  * Writes useful information concerning the matrix to a file. Should be
518  * executed after the matrix is factored.
519  *
520  * >>> Returns:
521  * One is returned if routine was successful, otherwise zero is returned.
522  * The calling function can query errno (the system global error variable)
523  * as to the reason why this routine failed.
524  *
525  * >>> Arguments:
526  * Matrix <input> (char *)
527  * Pointer to matrix.
528  * File <input> (char *)
529  * Name of file into which matrix is to be written.
530  * Label <input> (char *)
531  * String that is transferred to file and is used as a label.
532  *
533  * >>> Local variables:
534  * Data (RealNumber)
535  * The value of the matrix element being output.
536  * LargestElement (RealNumber)
537  * The largest element in the matrix.
538  * NumberOfElements (int)
539  * Number of nonzero elements in the matrix.
540  * pElement (ElementPtr)
541  * Pointer to an element in the matrix.
542  * pStatsFile (FILE *)
543  * File pointer to the statistics file.
544  * Size (int)
545  * The size of the matrix.
546  * SmallestElement (RealNumber)
547  * The smallest element in the matrix excluding zero elements.
548  */
549 
550 int spFileStats(char* eMatrix, char* File, char* Label)
551 {
552  MatrixPtr Matrix = (MatrixPtr)eMatrix;
553  int Size, I;
554  ElementPtr pElement;
555  int NumberOfElements;
556  RealNumber Data, LargestElement, SmallestElement;
557  FILE* pStatsFile;
558 
559  /* Begin `spFileStats'. */
561 
562  /* Open File in append mode. */
563  if ((pStatsFile = fopen(File, "a")) == NULL)
564  return 0;
565 
566  /* Output statistics. */
567  Size = Matrix->Size;
568  if (NOT Matrix->Factored)
569  fprintf(pStatsFile, "Matrix has not been factored.\n");
570  fprintf(pStatsFile, "||| Starting new matrix |||\n");
571  fprintf(pStatsFile, "%s\n", Label);
572  fprintf(pStatsFile, "Matrix is real.\n");
573  fprintf(pStatsFile, " Size = %d\n", Size);
574 
575  /* Search matrix. */
576  NumberOfElements = 0;
577  LargestElement = 0.0;
578  SmallestElement = DBL_MAX;
579 
580  for (I = 1; I <= Size; I++) {
581  pElement = Matrix->FirstInCol[I];
582  while (pElement != NULL) {
583  NumberOfElements++;
584  Data = ELEMENT_MAG(pElement);
585  if (Data > LargestElement)
586  LargestElement = Data;
587  if (Data < SmallestElement AND Data != 0.0)
588  SmallestElement = Data;
589  pElement = pElement->NextInCol;
590  }
591  }
592 
593  SmallestElement = MIN(SmallestElement, LargestElement);
594 
595  /* Output remaining statistics. */
596  fprintf(pStatsFile, " Initial number of elements = %d\n",
597  NumberOfElements - Matrix->Fillins);
598  fprintf(pStatsFile,
599  " Initial average number of elements per row = %lf\n",
600  (double)(NumberOfElements - Matrix->Fillins) / (double)Size);
601  fprintf(pStatsFile, " Fill-ins = %d\n", Matrix->Fillins);
602  fprintf(pStatsFile, " Average number of fill-ins per row = %lf%%\n",
603  (double)Matrix->Fillins / (double)Size);
604  fprintf(pStatsFile, " Total number of elements = %d\n",
605  NumberOfElements);
606  fprintf(pStatsFile, " Average number of elements per row = %lf\n",
607  (double)NumberOfElements / (double)Size);
608  fprintf(pStatsFile, " Density = %lf%%\n",
609  (double)(100.0 * NumberOfElements) / (double)(Size * Size));
610  fprintf(pStatsFile, " Relative Threshold = %e\n", Matrix->RelThreshold);
611  fprintf(pStatsFile, " Absolute Threshold = %e\n", Matrix->AbsThreshold);
612  fprintf(pStatsFile, " Largest Element = %e\n", LargestElement);
613  fprintf(pStatsFile, " Smallest Element = %e\n\n\n", SmallestElement);
614 
615  /* Close file. */
616  (void)fclose(pStatsFile);
617  return 1;
618 }
619 #endif /* DOCUMENTATION */
#define Label
Definition: _defines.h:157
#define MIN(a, b)
Definition: grids.h:32
printf
Definition: extdef.h:5
#define RHS(i)
Definition: multisplit.cpp:57
#define NULL
Definition: spdefs.h:105
#define FREE(ptr)
Definition: spdefs.h:177
#define OR
Definition: spdefs.h:101
struct MatrixFrame * MatrixPtr
Definition: spdefs.h:549
spREAL RealNumber
Definition: spdefs.h:202
#define CALLOC(ptr, type, number)
Definition: spdefs.h:187
spREAL * RealVector
Definition: spdefs.h:202
#define ASSERT(condition)
Definition: spdefs.h:151
#define AND
Definition: spdefs.h:100
#define NOT
Definition: spdefs.h:99
#define ELEMENT_MAG(ptr)
Definition: spdefs.h:134
#define IS_SPARSE(matrix)
Definition: spdefs.h:109
#define spNO_MEMORY
Definition: spmatrix.h:90
int spFileVector(char *eMatrix, char *File, RealVector RHS, RealVector iRHS)
Definition: spoutput.cpp:479
int spFileMatrix(char *eMatrix, char *File, char *Label, int Reordered, int Data, int Header)
Definition: spoutput.cpp:366
void spPrint(char *eMatrix, int PrintReordered, int Data, int Header)
Definition: spoutput.cpp:126
int spFileStats(char *eMatrix, char *File, char *Label)
Definition: spoutput.cpp:550
RealNumber Real
Definition: spdefs.h:245
struct MatrixElement * NextInCol
Definition: spdefs.h:249