projects / pspp / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
pivot-table: New functions for setting captions, etc.
[pspp] / src / language / stats / factor.c
diff --git a/src/language/stats/factor.c b/src/language/stats/factor.c
index a9c948e6369671ec6aa40f9a8103f91dbebb5ea4..3c1ce6094e5a2854be019066efdffae0e6afa618 100644 (file)
--- a/src/language/stats/factor.c
+++ b/src/language/stats/factor.c
@@ -1,5 +1,6 @@
 /* PSPP - a program for statistical analysis.
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2009, 2010, 2011 Free Software Foundation, Inc.
+   Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015,
+   2016, 2017 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -19,10 +20,12 @@
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
-#include <gsl/gsl_eigen.h> 
-#include <gsl/gsl_blas.h> 
+#include <gsl/gsl_eigen.h>
+#include <gsl/gsl_blas.h>
 #include <gsl/gsl_sort_vector.h>
 #include <gsl/gsl_sort_vector.h>
+#include <gsl/gsl_cdf.h>
 
 
+#include "data/any-reader.h"
 #include "data/casegrouper.h"
 #include "data/casereader.h"
 #include "data/casewriter.h"
 #include "data/casegrouper.h"
 #include "data/casereader.h"
 #include "data/casewriter.h"
@@ -34,14 +37,17 @@
 #include "language/lexer/lexer.h"
 #include "language/lexer/value-parser.h"
 #include "language/lexer/variable-parser.h"
 #include "language/lexer/lexer.h"
 #include "language/lexer/value-parser.h"
 #include "language/lexer/variable-parser.h"
+#include "language/data-io/file-handle.h"
+#include "language/data-io/matrix-reader.h"
+#include "libpspp/cast.h"
 #include "libpspp/message.h"
 #include "libpspp/misc.h"
 #include "math/correlation.h"
 #include "math/covariance.h"
 #include "math/moments.h"
 #include "libpspp/message.h"
 #include "libpspp/misc.h"
 #include "math/correlation.h"
 #include "math/covariance.h"
 #include "math/moments.h"
-#include "output/chart-item.h"
 #include "output/charts/scree.h"
 #include "output/charts/scree.h"
-#include "output/tab.h"
+#include "output/pivot-table.h"
+
 
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
 
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
@@ -85,7 +91,7 @@ enum print_opts
     PRINT_EXTRACTION  = 0x0100,
     PRINT_INITIAL     = 0x0200,
     PRINT_KMO         = 0x0400,
     PRINT_EXTRACTION  = 0x0100,
     PRINT_INITIAL     = 0x0200,
     PRINT_KMO         = 0x0400,
-    PRINT_REPR        = 0x0800, 
+    PRINT_REPR        = 0x0800,
     PRINT_FSCORE      = 0x1000
   };
 
     PRINT_FSCORE      = 0x1000
   };
 
@@ -94,18 +100,19 @@ enum rotation_type
     ROT_VARIMAX = 0,
     ROT_EQUAMAX,
     ROT_QUARTIMAX,
     ROT_VARIMAX = 0,
     ROT_EQUAMAX,
     ROT_QUARTIMAX,
+    ROT_PROMAX,
     ROT_NONE
   };
 
 typedef void (*rotation_coefficients) (double *x, double *y,
                                    double a, double b, double c, double d,
     ROT_NONE
   };
 
 typedef void (*rotation_coefficients) (double *x, double *y,
                                    double a, double b, double c, double d,
-                                   const gsl_matrix *loadings );
+                                   const gsl_matrix *loadings);
 
 
 static void
 varimax_coefficients (double *x, double *y,
                      double a, double b, double c, double d,
 
 
 static void
 varimax_coefficients (double *x, double *y,
                      double a, double b, double c, double d,
-                     const gsl_matrix *loadings )
+                     const gsl_matrix *loadings)
 {
   *x = d - 2 * a * b / loadings->size1;
   *y = c - (a * a - b * b) / loadings->size1;
 {
   *x = d - 2 * a * b / loadings->size1;
   *y = c - (a * a - b * b) / loadings->size1;
@@ -114,7 +121,7 @@ varimax_coefficients (double *x, double *y,
 static void
 equamax_coefficients (double *x, double *y,
                      double a, double b, double c, double d,
 static void
 equamax_coefficients (double *x, double *y,
                      double a, double b, double c, double d,
-                     const gsl_matrix *loadings )
+                     const gsl_matrix *loadings)
 {
   *x = d - loadings->size2 * a * b / loadings->size1;
   *y = c - loadings->size2 * (a * a - b * b) / (2 * loadings->size1);
 {
   *x = d - loadings->size2 * a * b / loadings->size1;
   *y = c - loadings->size2 * (a * a - b * b) / (2 * loadings->size1);
@@ -129,14 +136,66 @@ quartimax_coefficients (double *x, double *y,
   *y = c ;
 }
 
   *y = c ;
 }
 
-static const rotation_coefficients rotation_coeff[3] = {
+static const rotation_coefficients rotation_coeff[] = {
   varimax_coefficients,
   equamax_coefficients,
   varimax_coefficients,
   equamax_coefficients,
-  quartimax_coefficients
+  quartimax_coefficients,
+  varimax_coefficients  /* PROMAX is identical to VARIMAX */
 };
 
 
 };
 
 
-struct cmd_factor 
+/* return diag (C'C) ^ {-0.5} */
+static gsl_matrix *
+diag_rcp_sqrt (const gsl_matrix *C)
+{
+  int j;
+  gsl_matrix *d =  gsl_matrix_calloc (C->size1, C->size2);
+  gsl_matrix *r =  gsl_matrix_calloc (C->size1, C->size2);
+
+  assert (C->size1 == C->size2);
+
+  gsl_linalg_matmult_mod (C,  GSL_LINALG_MOD_TRANSPOSE,
+                         C,  GSL_LINALG_MOD_NONE,
+                         d);
+
+  for (j = 0 ; j < d->size2; ++j)
+    {
+      double e = gsl_matrix_get (d, j, j);
+      e = 1.0 / sqrt (e);
+      gsl_matrix_set (r, j, j, e);
+    }
+
+  gsl_matrix_free (d);
+
+  return r;
+}
+
+
+
+/* return diag ((C'C)^-1) ^ {-0.5} */
+static gsl_matrix *
+diag_rcp_inv_sqrt (const gsl_matrix *CCinv)
+{
+  int j;
+  gsl_matrix *r =  gsl_matrix_calloc (CCinv->size1, CCinv->size2);
+
+  assert (CCinv->size1 == CCinv->size2);
+
+  for (j = 0 ; j < CCinv->size2; ++j)
+    {
+      double e = gsl_matrix_get (CCinv, j, j);
+      e = 1.0 / sqrt (e);
+      gsl_matrix_set (r, j, j, e);
+    }
+
+  return r;
+}
+
+
+
+
+
+struct cmd_factor
 {
   size_t n_vars;
   const struct variable **vars;
 {
   size_t n_vars;
   const struct variable **vars;
@@ -150,12 +209,14 @@ struct cmd_factor
   enum extraction_method extraction;
   enum plot_opts plot;
   enum rotation_type rotation;
   enum extraction_method extraction;
   enum plot_opts plot;
   enum rotation_type rotation;
+  int rotation_iterations;
+  int promax_power;
 
   /* Extraction Criteria */
   int n_factors;
   double min_eigen;
   double econverge;
 
   /* Extraction Criteria */
   int n_factors;
   double min_eigen;
   double econverge;
-  int iterations;
+  int extraction_iterations;
 
   double rconverge;
 
 
   double rconverge;
 
@@ -164,13 +225,13 @@ struct cmd_factor
   bool sort;
 };
 
   bool sort;
 };
 
+
 struct idata
 {
   /* Intermediate values used in calculation */
 struct idata
 {
   /* Intermediate values used in calculation */
+  struct matrix_material mm;
 
 
-  const gsl_matrix *corr ;  /* The correlation matrix */
-  gsl_matrix *cov ;         /* The covariance matrix */
-  const gsl_matrix *n ;     /* Matrix of number of samples */
+  gsl_matrix *analysis_matrix; /* A pointer to either mm.corr or mm.cov */
 
   gsl_vector *eval ;  /* The eigenvalues */
   gsl_matrix *evec ;  /* The eigenvectors */
 
   gsl_vector *eval ;  /* The eigenvalues */
   gsl_matrix *evec ;  /* The eigenvectors */
@@ -178,12 +239,18 @@ struct idata
   int n_extractions;
 
   gsl_vector *msr ;  /* Multiple Squared Regressions */
   int n_extractions;
 
   gsl_vector *msr ;  /* Multiple Squared Regressions */
+
+  double detR;  /* The determinant of the correlation matrix */
+
+  gsl_matrix *ai_cov; /* The anti-image covariance matrix */
+  gsl_matrix *ai_cor; /* The anti-image correlation matrix */
+  struct covariance *cvm;
 };
 
 static struct idata *
 idata_alloc (size_t n_vars)
 {
 };
 
 static struct idata *
 idata_alloc (size_t n_vars)
 {
-  struct idata *id = xzalloc (sizeof (*id));
+  struct idata *id = XZALLOC (struct idata);
 
   id->n_extractions = 0;
   id->msr = gsl_vector_alloc (n_vars);
 
   id->n_extractions = 0;
   id->msr = gsl_vector_alloc (n_vars);
@@ -200,12 +267,106 @@ idata_free (struct idata *id)
   gsl_vector_free (id->msr);
   gsl_vector_free (id->eval);
   gsl_matrix_free (id->evec);
   gsl_vector_free (id->msr);
   gsl_vector_free (id->eval);
   gsl_matrix_free (id->evec);
-  if (id->cov != NULL)
-    gsl_matrix_free (id->cov);
+  gsl_matrix_free (id->ai_cov);
+  gsl_matrix_free (id->ai_cor);
 
   free (id);
 }
 
 
   free (id);
 }
 
+/* Return the sum of squares of all the elements in row J excluding column J */
+static double
+ssq_row_od_n (const gsl_matrix *m, int j)
+{
+  int i;
+  double ss = 0;
+  assert (m->size1 == m->size2);
+
+  assert (j < m->size1);
+
+  for (i = 0; i < m->size1; ++i)
+    {
+      if (i == j) continue;
+      ss += pow2 (gsl_matrix_get (m, i, j));
+    }
+
+  return ss;
+}
+
+/* Return the sum of squares of all the elements excluding row N */
+static double
+ssq_od_n (const gsl_matrix *m, int n)
+{
+  int i, j;
+  double ss = 0;
+  assert (m->size1 == m->size2);
+
+  assert (n < m->size1);
+
+  for (i = 0; i < m->size1; ++i)
+    {
+      for (j = 0; j < m->size2; ++j)
+       {
+         if (i == j) continue;
+         ss += pow2 (gsl_matrix_get (m, i, j));
+       }
+    }
+
+  return ss;
+}
+
+
+static gsl_matrix *
+anti_image_corr (const gsl_matrix *m, const struct idata *idata)
+{
+  int i, j;
+  gsl_matrix *a;
+  assert (m->size1 == m->size2);
+
+  a = gsl_matrix_alloc (m->size1, m->size2);
+
+  for (i = 0; i < m->size1; ++i)
+    {
+      for (j = 0; j < m->size2; ++j)
+        {
+          double *p = gsl_matrix_ptr (a, i, j);
+          *p = gsl_matrix_get (m, i, j);
+          *p /= sqrt (gsl_matrix_get (m, i, i) *
+                      gsl_matrix_get (m, j, j));
+        }
+    }
+
+  for (i = 0; i < m->size1; ++i)
+    {
+      double r = ssq_row_od_n (idata->mm.corr, i);
+      double u = ssq_row_od_n (a, i);
+      gsl_matrix_set (a, i, i, r / (r + u));
+    }
+
+  return a;
+}
+
+static gsl_matrix *
+anti_image_cov (const gsl_matrix *m)
+{
+  int i, j;
+  gsl_matrix *a;
+  assert (m->size1 == m->size2);
+
+  a = gsl_matrix_alloc (m->size1, m->size2);
+
+  for (i = 0; i < m->size1; ++i)
+    {
+      for (j = 0; j < m->size2; ++j)
+       {
+         double *p = gsl_matrix_ptr (a, i, j);
+         *p = gsl_matrix_get (m, i, j);
+         *p /= gsl_matrix_get (m, i, i);
+         *p /= gsl_matrix_get (m, j, j);
+       }
+    }
+
+  return a;
+}
 
 #if 0
 static void
 
 #if 0
 static void
@@ -221,7 +382,6 @@ dump_matrix (const gsl_matrix *m)
     }
 }
 
     }
 }
 
-
 static void
 dump_matrix_permute (const gsl_matrix *m, const gsl_permutation *p)
 {
 static void
 dump_matrix_permute (const gsl_matrix *m, const gsl_permutation *p)
 {
@@ -249,13 +409,13 @@ dump_vector (const gsl_vector *v)
 #endif
 
 
 #endif
 
 
-static int 
+static int
 n_extracted_factors (const struct cmd_factor *factor, struct idata *idata)
 {
   int i;
 n_extracted_factors (const struct cmd_factor *factor, struct idata *idata)
 {
   int i;
-  
+
   /* If there is a cached value, then return that. */
   /* If there is a cached value, then return that. */
-  if ( idata->n_extractions != 0)
+  if (idata->n_extractions != 0)
     return idata->n_extractions;
 
   /* Otherwise, if the number of factors has been explicitly requested,
     return idata->n_extractions;
 
   /* Otherwise, if the number of factors has been explicitly requested,
@@ -265,7 +425,7 @@ n_extracted_factors (const struct cmd_factor *factor, struct idata *idata)
       idata->n_extractions = factor->n_factors;
       goto finish;
     }
       idata->n_extractions = factor->n_factors;
       goto finish;
     }
-  
+
   /* Use the MIN_EIGEN setting. */
   for (i = 0 ; i < idata->eval->size; ++i)
     {
   /* Use the MIN_EIGEN setting. */
   for (i = 0 ; i < idata->eval->size; ++i)
     {
@@ -300,7 +460,7 @@ struct smr_workspace
 {
   /* Copy of the subject */
   gsl_matrix *m;
 {
   /* Copy of the subject */
   gsl_matrix *m;
-  
+
   gsl_matrix *inverse;
 
   gsl_permutation *perm;
   gsl_matrix *inverse;
 
   gsl_permutation *perm;
@@ -313,7 +473,7 @@ struct smr_workspace
 static struct smr_workspace *ws_create (const gsl_matrix *input)
 {
   struct smr_workspace *ws = xmalloc (sizeof (*ws));
 static struct smr_workspace *ws_create (const gsl_matrix *input)
 {
   struct smr_workspace *ws = xmalloc (sizeof (*ws));
-  
+
   ws->m = gsl_matrix_alloc (input->size1, input->size2);
   ws->inverse = gsl_matrix_calloc (input->size1 - 1, input->size2 - 1);
   ws->perm = gsl_permutation_alloc (input->size1 - 1);
   ws->m = gsl_matrix_alloc (input->size1, input->size2);
   ws->inverse = gsl_matrix_calloc (input->size1 - 1, input->size2 - 1);
   ws->perm = gsl_permutation_alloc (input->size1 - 1);
@@ -336,13 +496,13 @@ ws_destroy (struct smr_workspace *ws)
 }
 
 
 }
 
 
-/* 
+/*
    Return the square of the regression coefficient for VAR regressed against all other variables.
  */
 static double
 squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_workspace *ws)
 {
    Return the square of the regression coefficient for VAR regressed against all other variables.
  */
 static double
 squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_workspace *ws)
 {
-  /* For an explanation of what this is doing, see 
+  /* For an explanation of what this is doing, see
      http://www.visualstatistics.net/Visual%20Statistics%20Multimedia/multiple_regression_analysis.htm
   */
 
      http://www.visualstatistics.net/Visual%20Statistics%20Multimedia/multiple_regression_analysis.htm
   */
 
@@ -354,7 +514,7 @@ squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_worksp
   gsl_matrix_swap_rows (ws->m, 0, var);
   gsl_matrix_swap_columns (ws->m, 0, var);
 
   gsl_matrix_swap_rows (ws->m, 0, var);
   gsl_matrix_swap_columns (ws->m, 0, var);
 
-  rxx = gsl_matrix_submatrix (ws->m, 1, 1, ws->m->size1 - 1, ws->m->size1 - 1); 
+  rxx = gsl_matrix_submatrix (ws->m, 1, 1, ws->m->size1 - 1, ws->m->size1 - 1);
 
   gsl_linalg_LU_decomp (&rxx.matrix, ws->perm, &signum);
 
 
   gsl_linalg_LU_decomp (&rxx.matrix, ws->perm, &signum);
 
@@ -403,7 +563,7 @@ factor_matrix_workspace_alloc (size_t n, size_t nf)
   ws->eval = gsl_vector_alloc (n);
   ws->evec = gsl_matrix_alloc (n, n);
   ws->r  = gsl_matrix_alloc (n, n);
   ws->eval = gsl_vector_alloc (n);
   ws->evec = gsl_matrix_alloc (n, n);
   ws->r  = gsl_matrix_alloc (n, n);
-  
+
   return ws;
 }
 
   return ws;
 }
 
@@ -438,10 +598,10 @@ perm_shift_apply (gsl_permutation *target, const gsl_permutation *p,
 }
 
 
 }
 
 
-/* 
+/*
    Indirectly sort the rows of matrix INPUT, storing the sort order in PERM.
    The sort criteria are as follows:
    Indirectly sort the rows of matrix INPUT, storing the sort order in PERM.
    The sort criteria are as follows:
-   
+
    Rows are sorted on the first column, until the absolute value of an
    element in a subsequent column  is greater than that of the first
    column.  Thereafter, rows will be sorted on the second column,
    Rows are sorted on the first column, until the absolute value of an
    element in a subsequent column  is greater than that of the first
    column.  Thereafter, rows will be sorted on the second column,
@@ -474,7 +634,7 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm)
        }
     }
 
        }
     }
 
-  while (column_n < m && row_n < n) 
+  while (column_n < m && row_n < n)
     {
       gsl_vector_const_view columni = gsl_matrix_const_column (mat, column_n);
       gsl_sort_vector_index (p, &columni.vector);
     {
       gsl_vector_const_view columni = gsl_matrix_const_column (mat, column_n);
       gsl_sort_vector_index (p, &columni.vector);
@@ -483,8 +643,8 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm)
        {
          gsl_vector_view row = gsl_matrix_row (mat, p->data[n - 1 - i]);
          size_t maxindex = gsl_vector_max_index (&row.vector);
        {
          gsl_vector_view row = gsl_matrix_row (mat, p->data[n - 1 - i]);
          size_t maxindex = gsl_vector_max_index (&row.vector);
-         
-         if ( maxindex > column_n )
+
+         if (maxindex > column_n)
            break;
 
          /* All subsequent elements of this row, are of no interest.
            break;
 
          /* All subsequent elements of this row, are of no interest.
@@ -501,11 +661,11 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm)
 
   gsl_permutation_free (p);
   gsl_matrix_free (mat);
 
   gsl_permutation_free (p);
   gsl_matrix_free (mat);
-  
-  assert ( 0 == gsl_permutation_valid (perm));
+
+  assert (0 == gsl_permutation_valid (perm));
 
   /* We want the biggest value to be first */
 
   /* We want the biggest value to be first */
-  gsl_permutation_reverse (perm);    
+  gsl_permutation_reverse (perm);
 }
 
 
 }
 
 
@@ -539,7 +699,7 @@ clone_matrix (const gsl_matrix *m)
 }
 
 
 }
 
 
-static double 
+static double
 initial_sv (const gsl_matrix *fm)
 {
   int j, k;
 initial_sv (const gsl_matrix *fm)
 {
   int j, k;
@@ -559,7 +719,7 @@ initial_sv (const gsl_matrix *fm)
          l4s += lambda_4;
          l2s += lambda_sq;
        }
          l4s += lambda_4;
          l2s += lambda_sq;
        }
-      sv += ( fm->size1 * l4s - (l2s * l2s) ) / (fm->size1 * fm->size1 );
+      sv += (fm->size1 * l4s - (l2s * l2s)) / (fm->size1 * fm->size1);
     }
   return sv;
 }
     }
   return sv;
 }
@@ -568,7 +728,9 @@ static void
 rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
        const gsl_vector *communalities,
        gsl_matrix *result,
 rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
        const gsl_vector *communalities,
        gsl_matrix *result,
-       gsl_vector *rotated_loadings
+       gsl_vector *rotated_loadings,
+       gsl_matrix *pattern_matrix,
+       gsl_matrix *factor_correlation_matrix
        )
 {
   int j, k;
        )
 {
   int j, k;
@@ -605,7 +767,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
   /* Now perform the rotation iterations */
 
   prev_sv = initial_sv (normalised);
   /* Now perform the rotation iterations */
 
   prev_sv = initial_sv (normalised);
-  for (i = 0 ; i < cf->iterations ; ++i)
+  for (i = 0 ; i < cf->rotation_iterations ; ++i)
     {
       double sv = 0.0;
       for (j = 0 ; j < normalised->size2; ++j)
     {
       double sv = 0.0;
       for (j = 0 ; j < normalised->size2; ++j)
@@ -628,7 +790,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
                {
                  double jv = gsl_matrix_get (normalised, p, j);
                  double kv = gsl_matrix_get (normalised, p, k);
                {
                  double jv = gsl_matrix_get (normalised, p, j);
                  double kv = gsl_matrix_get (normalised, p, k);
-             
+
                  double u = jv * jv - kv * kv;
                  double v = 2 * jv * kv;
                  a += u;
                  double u = jv * jv - kv * kv;
                  double v = 2 * jv * kv;
                  a += u;
@@ -642,7 +804,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
              phi = atan2 (x,  y) / 4.0 ;
 
              /* Don't bother rotating if the angle is small */
              phi = atan2 (x,  y) / 4.0 ;
 
              /* Don't bother rotating if the angle is small */
-             if ( fabs (sin (phi) ) <= pow (10.0, -15.0))
+             if (fabs (sin (phi)) <= pow (10.0, -15.0))
                  continue;
 
              for (p = 0; p < normalised->size1; ++p)
                  continue;
 
              for (p = 0; p < normalised->size1; ++p)
@@ -662,10 +824,10 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
                l2s += lambda_sq;
              }
            }
                l2s += lambda_sq;
              }
            }
-         sv += ( normalised->size1 * l4s - (l2s * l2s) ) / (normalised->size1 * normalised->size1 );
+         sv += (normalised->size1 * l4s - (l2s * l2s)) / (normalised->size1 * normalised->size1);
        }
 
        }
 
-      if ( fabs (sv - prev_sv) <= cf->rconverge)
+      if (fabs (sv - prev_sv) <= cf->rconverge)
        break;
 
       prev_sv = sv;
        break;
 
       prev_sv = sv;
@@ -675,6 +837,143 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
                  h_sqrt, normalised,  0.0,   result);
 
   gsl_matrix_free (h_sqrt);
                  h_sqrt, normalised,  0.0,   result);
 
   gsl_matrix_free (h_sqrt);
+  gsl_matrix_free (normalised);
+
+  if (cf->rotation == ROT_PROMAX)
+    {
+      /* general purpose m by m matrix, where m is the number of factors */
+      gsl_matrix *mm1 =  gsl_matrix_calloc (unrot->size2, unrot->size2);
+      gsl_matrix *mm2 =  gsl_matrix_calloc (unrot->size2, unrot->size2);
+
+      /* general purpose m by p matrix, where p is the number of variables */
+      gsl_matrix *mp1 =  gsl_matrix_calloc (unrot->size2, unrot->size1);
+
+      gsl_matrix *pm1 =  gsl_matrix_calloc (unrot->size1, unrot->size2);
+
+      gsl_permutation *perm = gsl_permutation_alloc (unrot->size2);
+
+      int signum;
+
+      int i, j;
+
+      /* The following variables follow the notation by SPSS Statistical Algorithms
+        page 342 */
+      gsl_matrix *L =  gsl_matrix_calloc (unrot->size2, unrot->size2);
+      gsl_matrix *P = clone_matrix (result);
+      gsl_matrix *D ;
+      gsl_matrix *Q ;
+
+
+      /* Vector of length p containing (indexed by i)
+        \Sum^m_j {\lambda^2_{ij}} */
+      gsl_vector *rssq = gsl_vector_calloc (unrot->size1);
+
+      for (i = 0; i < P->size1; ++i)
+       {
+         double sum = 0;
+         for (j = 0; j < P->size2; ++j)
+           {
+             sum += gsl_matrix_get (result, i, j)
+               * gsl_matrix_get (result, i, j);
+
+           }
+
+         gsl_vector_set (rssq, i, sqrt (sum));
+       }
+
+      for (i = 0; i < P->size1; ++i)
+       {
+         for (j = 0; j < P->size2; ++j)
+           {
+             double l = gsl_matrix_get (result, i, j);
+             double r = gsl_vector_get (rssq, i);
+             gsl_matrix_set (P, i, j, pow (fabs (l / r), cf->promax_power + 1) * r / l);
+           }
+       }
+
+      gsl_vector_free (rssq);
+
+      gsl_linalg_matmult_mod (result,
+                             GSL_LINALG_MOD_TRANSPOSE,
+                             result,
+                             GSL_LINALG_MOD_NONE,
+                             mm1);
+
+      gsl_linalg_LU_decomp (mm1, perm, &signum);
+      gsl_linalg_LU_invert (mm1, perm, mm2);
+
+      gsl_linalg_matmult_mod (mm2,   GSL_LINALG_MOD_NONE,
+                             result,  GSL_LINALG_MOD_TRANSPOSE,
+                             mp1);
+
+      gsl_linalg_matmult_mod (mp1, GSL_LINALG_MOD_NONE,
+                             P,   GSL_LINALG_MOD_NONE,
+                             L);
+
+      D = diag_rcp_sqrt (L);
+      Q = gsl_matrix_calloc (unrot->size2, unrot->size2);
+
+      gsl_linalg_matmult_mod (L, GSL_LINALG_MOD_NONE,
+                             D, GSL_LINALG_MOD_NONE,
+                             Q);
+
+      gsl_matrix *QQinv = gsl_matrix_calloc (unrot->size2, unrot->size2);
+
+      gsl_linalg_matmult_mod (Q, GSL_LINALG_MOD_TRANSPOSE,
+                             Q,  GSL_LINALG_MOD_NONE,
+                             QQinv);
+
+      gsl_linalg_cholesky_decomp (QQinv);
+      gsl_linalg_cholesky_invert (QQinv);
+
+
+      gsl_matrix *C = diag_rcp_inv_sqrt (QQinv);
+      gsl_matrix *Cinv =  clone_matrix (C);
+
+      gsl_linalg_cholesky_decomp (Cinv);
+      gsl_linalg_cholesky_invert (Cinv);
+
+
+      gsl_linalg_matmult_mod (result, GSL_LINALG_MOD_NONE,
+                             Q,      GSL_LINALG_MOD_NONE,
+                             pm1);
+
+      gsl_linalg_matmult_mod (pm1,      GSL_LINALG_MOD_NONE,
+                             Cinv,         GSL_LINALG_MOD_NONE,
+                             pattern_matrix);
+
+
+      gsl_linalg_matmult_mod (C,      GSL_LINALG_MOD_NONE,
+                             QQinv,  GSL_LINALG_MOD_NONE,
+                             mm1);
+
+      gsl_linalg_matmult_mod (mm1,      GSL_LINALG_MOD_NONE,
+                             C,  GSL_LINALG_MOD_TRANSPOSE,
+                             factor_correlation_matrix);
+
+      gsl_linalg_matmult_mod (pattern_matrix,      GSL_LINALG_MOD_NONE,
+                             factor_correlation_matrix,  GSL_LINALG_MOD_NONE,
+                             pm1);
+
+      gsl_matrix_memcpy (result, pm1);
+
+
+      gsl_matrix_free (QQinv);
+      gsl_matrix_free (C);
+      gsl_matrix_free (Cinv);
+
+      gsl_matrix_free (D);
+      gsl_matrix_free (Q);
+      gsl_matrix_free (L);
+      gsl_matrix_free (P);
+
+      gsl_permutation_free (perm);
+
+      gsl_matrix_free (mm1);
+      gsl_matrix_free (mm2);
+      gsl_matrix_free (mp1);
+      gsl_matrix_free (pm1);
+    }
 
 
   /* reflect negative sums and populate the rotated loadings vector*/
 
 
   /* reflect negative sums and populate the rotated loadings vector*/
@@ -686,12 +985,12 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
        {
          double s = gsl_matrix_get (result, j, i);
          ssq += s * s;
        {
          double s = gsl_matrix_get (result, j, i);
          ssq += s * s;
-         sum += gsl_matrix_get (result, j, i);
+         sum += s;
        }
 
       gsl_vector_set (rotated_loadings, i, ssq);
 
        }
 
       gsl_vector_set (rotated_loadings, i, ssq);
 
-      if ( sum < 0 )
+      if (sum < 0)
        for (j = 0 ; j < result->size1; ++j)
          {
            double *lambda = gsl_matrix_ptr (result, j, i);
        for (j = 0 ; j < result->size1; ++j)
          {
            double *lambda = gsl_matrix_ptr (result, j, i);
@@ -707,7 +1006,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
   WS is a pointer to a structure which must have been initialised with factor_matrix_workspace_init.
  */
 static void
   WS is a pointer to a structure which must have been initialised with factor_matrix_workspace_init.
  */
 static void
-iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors, 
+iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors,
                       struct factor_matrix_workspace *ws)
 {
   size_t i;
                       struct factor_matrix_workspace *ws)
 {
   size_t i;
@@ -755,13 +1054,15 @@ iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matri
 
 static bool run_factor (struct dataset *ds, const struct cmd_factor *factor);
 
 
 static bool run_factor (struct dataset *ds, const struct cmd_factor *factor);
 
+static void do_factor_by_matrix (const struct cmd_factor *factor, struct idata *idata);
+
+
 
 int
 cmd_factor (struct lexer *lexer, struct dataset *ds)
 {
 
 int
 cmd_factor (struct lexer *lexer, struct dataset *ds)
 {
-  bool extraction_seen = false;
-  const struct dictionary *dict = dataset_dict (ds);
-
+  struct dictionary *dict = NULL;
+  int n_iterations = 25;
   struct cmd_factor factor;
   factor.n_vars = 0;
   factor.vars = NULL;
   struct cmd_factor factor;
   factor.n_vars = 0;
   factor.vars = NULL;
@@ -772,39 +1073,120 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
   factor.extraction = EXTRACTION_PC;
   factor.n_factors = 0;
   factor.min_eigen = SYSMIS;
   factor.extraction = EXTRACTION_PC;
   factor.n_factors = 0;
   factor.min_eigen = SYSMIS;
-  factor.iterations = 25;
+  factor.extraction_iterations = 25;
+  factor.rotation_iterations = 25;
   factor.econverge = 0.001;
 
   factor.blank = 0;
   factor.sort = false;
   factor.plot = 0;
   factor.rotation = ROT_VARIMAX;
   factor.econverge = 0.001;
 
   factor.blank = 0;
   factor.sort = false;
   factor.plot = 0;
   factor.rotation = ROT_VARIMAX;
+  factor.wv = NULL;
 
   factor.rconverge = 0.0001;
 
 
   factor.rconverge = 0.0001;
 
-  factor.wv = dict_get_weight (dict);
-
   lex_match (lexer, T_SLASH);
 
   lex_match (lexer, T_SLASH);
 
-  if (!lex_force_match_id (lexer, "VARIABLES"))
+  struct matrix_reader *mr = NULL;
+  struct casereader *matrix_reader = NULL;
+
+  if (lex_match_id (lexer, "VARIABLES"))
     {
     {
-      goto error;
+      lex_match (lexer, T_EQUALS);
+      dict = dataset_dict (ds);
+      factor.wv = dict_get_weight (dict);
+
+      if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars,
+                                 PV_NO_DUPLICATE | PV_NUMERIC))
+       goto error;
     }
     }
+  else if (lex_match_id (lexer, "MATRIX"))
+    {
+      lex_match (lexer, T_EQUALS);
+      if (! lex_force_match_id (lexer, "IN"))
+       goto error;
+      if (!lex_force_match (lexer, T_LPAREN))
+       {
+         goto error;
+       }
+      if (lex_match_id (lexer, "CORR"))
+       {
+       }
+      else if (lex_match_id (lexer, "COV"))
+       {
+       }
+      else
+       {
+         lex_error (lexer, _("Matrix input for %s must be either COV or CORR"), "FACTOR");
+         goto error;
+       }
+      if (! lex_force_match (lexer, T_EQUALS))
+       goto error;
+      if (lex_match (lexer, T_ASTERISK))
+       {
+         dict = dataset_dict (ds);
+         matrix_reader = casereader_clone (dataset_source (ds));
+       }
+      else
+       {
+         struct file_handle *fh = fh_parse (lexer, FH_REF_FILE, NULL);
+         if (fh == NULL)
+           goto error;
 
 
-  lex_match (lexer, T_EQUALS);
+         matrix_reader
+           = any_reader_open_and_decode (fh, NULL, &dict, NULL);
 
 
-  if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars,
-                             PV_NO_DUPLICATE | PV_NUMERIC))
-    goto error;
+         if (! (matrix_reader && dict))
+           {
+             goto error;
+           }
+       }
 
 
-  if (factor.n_vars < 2)
-    msg (MW, _("Factor analysis on a single variable is not useful."));
+      if (! lex_force_match (lexer, T_RPAREN))
+       goto error;
+
+      mr = matrix_reader_create (dict, matrix_reader);
+      factor.vars = xmemdup (mr->cvars, mr->n_cvars * sizeof *mr->cvars);
+      factor.n_vars = mr->n_cvars;
+    }
+  else
+    {
+      goto error;
+    }
 
   while (lex_token (lexer) != T_ENDCMD)
     {
       lex_match (lexer, T_SLASH);
 
 
   while (lex_token (lexer) != T_ENDCMD)
     {
       lex_match (lexer, T_SLASH);
 
-      if (lex_match_id (lexer, "PLOT"))
+      if (lex_match_id (lexer, "ANALYSIS"))
+        {
+          struct const_var_set *vs;
+          const struct variable **vars;
+          size_t n_vars;
+          bool ok;
+
+          lex_match (lexer, T_EQUALS);
+
+          vs = const_var_set_create_from_array (factor.vars, factor.n_vars);
+          ok = parse_const_var_set_vars (lexer, vs, &vars, &n_vars,
+                                         PV_NO_DUPLICATE | PV_NUMERIC);
+          const_var_set_destroy (vs);
+
+          if (!ok)
+            goto error;
+
+          free (factor.vars);
+          factor.vars = vars;
+          factor.n_vars = n_vars;
+
+          if (mr)
+            {
+              free (mr->cvars);
+              mr->cvars = xmemdup (vars, n_vars * sizeof *vars);
+              mr->n_cvars = n_vars;
+            }
+        }
+      else if (lex_match_id (lexer, "PLOT"))
        {
           lex_match (lexer, T_EQUALS);
           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
        {
           lex_match (lexer, T_EQUALS);
           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
@@ -863,6 +1245,19 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                {
                  factor.rotation = ROT_QUARTIMAX;
                }
                {
                  factor.rotation = ROT_QUARTIMAX;
                }
+             else if (lex_match_id (lexer, "PROMAX"))
+               {
+                 factor.promax_power = 5;
+                 if (lex_match (lexer, T_LPAREN)
+                      && lex_force_int (lexer))
+                   {
+                     factor.promax_power = lex_integer (lexer);
+                     lex_get (lexer);
+                     if (! lex_force_match (lexer, T_RPAREN))
+                       goto error;
+                   }
+                 factor.rotation = ROT_PROMAX;
+               }
              else if (lex_match_id (lexer, "NOROTATE"))
                {
                  factor.rotation = ROT_NONE;
              else if (lex_match_id (lexer, "NOROTATE"))
                {
                  factor.rotation = ROT_NONE;
@@ -873,6 +1268,7 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                  goto error;
                }
            }
                  goto error;
                }
            }
+          factor.rotation_iterations = n_iterations;
        }
       else if (lex_match_id (lexer, "CRITERIA"))
        {
        }
       else if (lex_match_id (lexer, "CRITERIA"))
        {
@@ -881,59 +1277,64 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
            {
              if (lex_match_id (lexer, "FACTORS"))
                {
            {
              if (lex_match_id (lexer, "FACTORS"))
                {
-                 if ( lex_force_match (lexer, T_LPAREN))
+                 if (lex_force_match (lexer, T_LPAREN)
+                       && lex_force_int (lexer))
                    {
                    {
-                     lex_force_int (lexer);
                      factor.n_factors = lex_integer (lexer);
                      lex_get (lexer);
                      factor.n_factors = lex_integer (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, T_RPAREN);
+                     if (! lex_force_match (lexer, T_RPAREN))
+                       goto error;
                    }
                }
              else if (lex_match_id (lexer, "MINEIGEN"))
                {
                    }
                }
              else if (lex_match_id (lexer, "MINEIGEN"))
                {
-                 if ( lex_force_match (lexer, T_LPAREN))
+                 if (lex_force_match (lexer, T_LPAREN)
+                       && lex_force_num (lexer))
                    {
                    {
-                     lex_force_num (lexer);
                      factor.min_eigen = lex_number (lexer);
                      lex_get (lexer);
                      factor.min_eigen = lex_number (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, T_RPAREN);
+                     if (! lex_force_match (lexer, T_RPAREN))
+                       goto error;
                    }
                }
              else if (lex_match_id (lexer, "ECONVERGE"))
                {
                    }
                }
              else if (lex_match_id (lexer, "ECONVERGE"))
                {
-                 if ( lex_force_match (lexer, T_LPAREN))
+                 if (lex_force_match (lexer, T_LPAREN)
+                       && lex_force_num (lexer))
                    {
                    {
-                     lex_force_num (lexer);
                      factor.econverge = lex_number (lexer);
                      lex_get (lexer);
                      factor.econverge = lex_number (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, T_RPAREN);
+                     if (! lex_force_match (lexer, T_RPAREN))
+                       goto error;
                    }
                }
              else if (lex_match_id (lexer, "RCONVERGE"))
                    }
                }
              else if (lex_match_id (lexer, "RCONVERGE"))
-               {
-                 if ( lex_force_match (lexer, T_LPAREN))
-                   {
-                     lex_force_num (lexer);
-                     factor.rconverge = lex_number (lexer);
-                     lex_get (lexer);
-                     lex_force_match (lexer, T_RPAREN);
-                   }
+                {
+                  if (lex_force_match (lexer, T_LPAREN)
+                      && lex_force_num (lexer))
+                    {
+                      factor.rconverge = lex_number (lexer);
+                      lex_get (lexer);
+                      if (! lex_force_match (lexer, T_RPAREN))
+                       goto error;
+                    }
                }
              else if (lex_match_id (lexer, "ITERATE"))
                {
                }
              else if (lex_match_id (lexer, "ITERATE"))
                {
-                 if ( lex_force_match (lexer, T_LPAREN))
+                 if (lex_force_match (lexer, T_LPAREN)
+                      && lex_force_int_range (lexer, "ITERATE", 0, INT_MAX))
                    {
                    {
-                     lex_force_int (lexer);
-                     factor.iterations = lex_integer (lexer);
+                     n_iterations = lex_integer (lexer);
                      lex_get (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, T_RPAREN);
+                     if (! lex_force_match (lexer, T_RPAREN))
+                       goto error;
                    }
                }
              else if (lex_match_id (lexer, "DEFAULT"))
                {
                  factor.n_factors = 0;
                  factor.min_eigen = 1;
                    }
                }
              else if (lex_match_id (lexer, "DEFAULT"))
                {
                  factor.n_factors = 0;
                  factor.min_eigen = 1;
-                 factor.iterations = 25;
+                 n_iterations = 25;
                }
              else
                {
                }
              else
                {
@@ -944,7 +1345,6 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
        }
       else if (lex_match_id (lexer, "EXTRACTION"))
        {
        }
       else if (lex_match_id (lexer, "EXTRACTION"))
        {
-         extraction_seen = true;
           lex_match (lexer, T_EQUALS);
           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
            {
           lex_match (lexer, T_EQUALS);
           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
            {
@@ -970,6 +1370,7 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                  goto error;
                }
            }
                  goto error;
                }
            }
+          factor.extraction_iterations = n_iterations;
        }
       else if (lex_match_id (lexer, "FORMAT"))
        {
        }
       else if (lex_match_id (lexer, "FORMAT"))
        {
@@ -982,12 +1383,13 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                }
              else if (lex_match_id (lexer, "BLANK"))
                {
                }
              else if (lex_match_id (lexer, "BLANK"))
                {
-                 if ( lex_force_match (lexer, T_LPAREN))
+                 if (lex_force_match (lexer, T_LPAREN)
+                       && lex_force_num (lexer))
                    {
                    {
-                     lex_force_num (lexer);
                      factor.blank = lex_number (lexer);
                      lex_get (lexer);
                      factor.blank = lex_number (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, T_RPAREN);
+                     if (! lex_force_match (lexer, T_RPAREN))
+                       goto error;
                    }
                }
              else if (lex_match_id (lexer, "DEFAULT"))
                    }
                }
              else if (lex_match_id (lexer, "DEFAULT"))
@@ -1020,10 +1422,11 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
              else if (lex_match_id (lexer, "INV"))
                {
                }
              else if (lex_match_id (lexer, "INV"))
                {
                }
+#endif
              else if (lex_match_id (lexer, "AIC"))
                {
              else if (lex_match_id (lexer, "AIC"))
                {
+                 factor.print |= PRINT_AIC;
                }
                }
-#endif
              else if (lex_match_id (lexer, "SIG"))
                {
                  factor.print |= PRINT_SIG;
              else if (lex_match_id (lexer, "SIG"))
                {
                  factor.print |= PRINT_SIG;
@@ -1032,11 +1435,10 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                {
                  factor.print |= PRINT_CORRELATION;
                }
                {
                  factor.print |= PRINT_CORRELATION;
                }
-#if FACTOR_FULLY_IMPLEMENTED
              else if (lex_match_id (lexer, "COVARIANCE"))
                {
              else if (lex_match_id (lexer, "COVARIANCE"))
                {
+                 factor.print |= PRINT_COVARIANCE;
                }
                }
-#endif
              else if (lex_match_id (lexer, "ROTATION"))
                {
                  factor.print |= PRINT_ROTATION;
              else if (lex_match_id (lexer, "ROTATION"))
                {
                  factor.print |= PRINT_ROTATION;
@@ -1049,10 +1451,11 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                {
                  factor.print |= PRINT_INITIAL;
                }
                {
                  factor.print |= PRINT_INITIAL;
                }
-#if FACTOR_FULLY_IMPLEMENTED
              else if (lex_match_id (lexer, "KMO"))
                {
              else if (lex_match_id (lexer, "KMO"))
                {
+                 factor.print |= PRINT_KMO;
                }
                }
+#if FACTOR_FULLY_IMPLEMENTED
              else if (lex_match_id (lexer, "REPR"))
                {
                }
              else if (lex_match_id (lexer, "REPR"))
                {
                }
@@ -1116,16 +1519,46 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
        }
     }
 
        }
     }
 
-  if ( factor.rotation == ROT_NONE )
+  if (factor.rotation == ROT_NONE)
     factor.print &= ~PRINT_ROTATION;
 
     factor.print &= ~PRINT_ROTATION;
 
-  if ( ! run_factor (ds, &factor)) 
-    goto error;
+  if (factor.n_vars < 2)
+    msg (MW, _("Factor analysis on a single variable is not useful."));
+
+  if (factor.n_vars < 1)
+    {
+      msg (ME, _("Factor analysis without variables is not possible."));
+      goto error;
+    }
 
 
+  if (matrix_reader)
+    {
+      struct idata *id = idata_alloc (factor.n_vars);
+
+      while (matrix_reader_next (&id->mm, mr, NULL))
+       {
+         do_factor_by_matrix (&factor, id);
+
+          gsl_matrix_free (id->ai_cov);
+          id->ai_cov = NULL;
+          gsl_matrix_free (id->ai_cor);
+          id->ai_cor = NULL;
+
+          matrix_material_uninit (&id->mm);
+       }
+
+      idata_free (id);
+    }
+  else
+    if (! run_factor (ds, &factor))
+      goto error;
+
+  matrix_reader_destroy (mr);
   free (factor.vars);
   return CMD_SUCCESS;
 
  error:
   free (factor.vars);
   return CMD_SUCCESS;
 
  error:
+  matrix_reader_destroy (mr);
   free (factor.vars);
   return CMD_FAILURE;
 }
   free (factor.vars);
   return CMD_FAILURE;
 }
@@ -1144,7 +1577,7 @@ run_factor (struct dataset *ds, const struct cmd_factor *factor)
 
   while (casegrouper_get_next_group (grouper, &group))
     {
 
   while (casegrouper_get_next_group (grouper, &group))
     {
-      if ( factor->missing_type == MISS_LISTWISE )
+      if (factor->missing_type == MISS_LISTWISE)
        group  = casereader_create_filter_missing (group, factor->vars, factor->n_vars,
                                                   factor->exclude,
                                                   NULL,  NULL);
        group  = casereader_create_filter_missing (group, factor->vars, factor->n_vars,
                                                   factor->exclude,
                                                   NULL,  NULL);
@@ -1185,19 +1618,19 @@ the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_fa
 
 /* Return the communality of variable N, calculated to N_FACTORS */
 static double
 
 /* Return the communality of variable N, calculated to N_FACTORS */
 static double
-communality (struct idata *idata, int n, int n_factors)
+communality (const struct idata *idata, int n, int n_factors)
 {
   return the_communality (idata->evec, idata->eval, n, n_factors);
 }
 
 
 static void
 {
   return the_communality (idata->evec, idata->eval, n, n_factors);
 }
 
 
 static void
-show_scree (const struct cmd_factor *f, struct idata *idata)
+show_scree (const struct cmd_factor *f, const struct idata *idata)
 {
   struct scree *s;
   const char *label ;
 
 {
   struct scree *s;
   const char *label ;
 
-  if ( !(f->plot & PLOT_SCREE) )
+  if (!(f->plot & PLOT_SCREE))
     return;
 
 
     return;
 
 
@@ -1212,583 +1645,535 @@ static void
 show_communalities (const struct cmd_factor * factor,
                    const gsl_vector *initial, const gsl_vector *extracted)
 {
 show_communalities (const struct cmd_factor * factor,
                    const gsl_vector *initial, const gsl_vector *extracted)
 {
-  int i;
-  int c = 0;
-  const int heading_columns = 1;
-  int nc = heading_columns;
-  const int heading_rows = 1;
-  const int nr = heading_rows + factor->n_vars;
-  struct tab_table *t;
-
-  if (factor->print & PRINT_EXTRACTION)
-    nc++;
-
-  if (factor->print & PRINT_INITIAL)
-    nc++;
-
-  /* No point having a table with only headings */
-  if (nc <= 1)
+  if (!(factor->print & (PRINT_INITIAL | PRINT_EXTRACTION)))
     return;
 
     return;
 
-  t = tab_create (nc, nr);
+  struct pivot_table *table = pivot_table_create (N_("Communalities"));
 
 
-  tab_title (t, _("Communalities"));
-
-  tab_headers (t, heading_columns, 0, heading_rows, 0);
-
-  c = 1;
+  struct pivot_dimension *communalities = pivot_dimension_create (
+    table, PIVOT_AXIS_COLUMN, N_("Communalities"));
   if (factor->print & PRINT_INITIAL)
   if (factor->print & PRINT_INITIAL)
-    tab_text (t, c++, 0, TAB_CENTER | TAT_TITLE, _("Initial"));
-
+    pivot_category_create_leaves (communalities->root, N_("Initial"));
   if (factor->print & PRINT_EXTRACTION)
   if (factor->print & PRINT_EXTRACTION)
-    tab_text (t, c++, 0, TAB_CENTER | TAT_TITLE, _("Extraction"));
-
-  /* Outline the box */
-  tab_box (t,
-          TAL_2, TAL_2,
-          -1, -1,
-          0, 0,
-          nc - 1, nr - 1);
-
-  /* Vertical lines */
-  tab_box (t,
-          -1, -1,
-          -1, TAL_1,
-          heading_columns, 0,
-          nc - 1, nr - 1);
-
-  tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
-  tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
-
-  for (i = 0 ; i < factor->n_vars; ++i)
+    pivot_category_create_leaves (communalities->root, N_("Extraction"));
+
+  struct pivot_dimension *variables = pivot_dimension_create (
+    table, PIVOT_AXIS_ROW, N_("Variables"));
+
+  for (size_t i = 0 ; i < factor->n_vars; ++i)
     {
     {
-      c = 0;
-      tab_text (t, c++, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[i]));
+      int row = pivot_category_create_leaf (
+        variables->root, pivot_value_new_variable (factor->vars[i]));
 
 
+      int col = 0;
       if (factor->print & PRINT_INITIAL)
       if (factor->print & PRINT_INITIAL)
-       tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL);
-
+        pivot_table_put2 (table, col++, row, pivot_value_new_number (
+                            gsl_vector_get (initial, i)));
       if (factor->print & PRINT_EXTRACTION)
       if (factor->print & PRINT_EXTRACTION)
-       tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL);
+        pivot_table_put2 (table, col++, row, pivot_value_new_number (
+                            gsl_vector_get (extracted, i)));
     }
 
     }
 
-  tab_submit (t);
+  pivot_table_submit (table);
 }
 
 }
 
+static struct pivot_dimension *
+create_numeric_dimension (struct pivot_table *table,
+                          enum pivot_axis_type axis_type, const char *name,
+                          size_t n, bool show_label)
+{
+  struct pivot_dimension *d = pivot_dimension_create (table, axis_type, name);
+  d->root->show_label = show_label;
+  for (int i = 0 ; i < n; ++i)
+    pivot_category_create_leaf (d->root, pivot_value_new_integer (i + 1));
+  return d;
+}
 
 static void
 
 static void
-show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const char *title, const gsl_matrix *fm)
+show_factor_matrix (const struct cmd_factor *factor, const struct idata *idata, const char *title, const gsl_matrix *fm)
 {
 {
-  int i;
-  const int n_factors = idata->n_extractions;
-
-  const int heading_columns = 1;
-  const int heading_rows = 2;
-  const int nr = heading_rows + factor->n_vars;
-  const int nc = heading_columns + n_factors;
-  gsl_permutation *perm;
-
-  struct tab_table *t = tab_create (nc, nr);
-
-  /* 
-  if ( factor->extraction == EXTRACTION_PC )
-    tab_title (t, _("Component Matrix"));
-  else 
-    tab_title (t, _("Factor Matrix"));
-  */
-
-  tab_title (t, "%s", title);
-
-  tab_headers (t, heading_columns, 0, heading_rows, 0);
-
-  if ( factor->extraction == EXTRACTION_PC )
-    tab_joint_text (t,
-                   1, 0,
-                   nc - 1, 0,
-                   TAB_CENTER | TAT_TITLE, _("Component"));
-  else
-    tab_joint_text (t,
-                   1, 0,
-                   nc - 1, 0,
-                   TAB_CENTER | TAT_TITLE, _("Factor"));
-
-
-  tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
-
-
-  /* Outline the box */
-  tab_box (t,
-          TAL_2, TAL_2,
-          -1, -1,
-          0, 0,
-          nc - 1, nr - 1);
+  struct pivot_table *table = pivot_table_create (title);
 
 
-  /* Vertical lines */
-  tab_box (t,
-          -1, -1,
-          -1, TAL_1,
-          heading_columns, 1,
-          nc - 1, nr - 1);
-
-  tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
-  tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
+  const int n_factors = idata->n_extractions;
+  create_numeric_dimension (
+    table, PIVOT_AXIS_COLUMN,
+    factor->extraction == EXTRACTION_PC ? N_("Component") : N_("Factor"),
+    n_factors, true);
 
 
+  struct pivot_dimension *variables = pivot_dimension_create (
+    table, PIVOT_AXIS_ROW, N_("Variables"));
 
   /* Initialise to the identity permutation */
 
   /* Initialise to the identity permutation */
-  perm = gsl_permutation_calloc (factor->n_vars);
+  gsl_permutation *perm = gsl_permutation_calloc (factor->n_vars);
 
 
-  if ( factor->sort)
+  if (factor->sort)
     sort_matrix_indirect (fm, perm);
 
     sort_matrix_indirect (fm, perm);
 
-  for (i = 0 ; i < n_factors; ++i)
+  for (size_t i = 0 ; i < factor->n_vars; ++i)
     {
     {
-      tab_text_format (t, heading_columns + i, 1, TAB_CENTER | TAT_TITLE, _("%d"), i + 1);
-    }
-
-  for (i = 0 ; i < factor->n_vars; ++i)
-    {
-      int j;
       const int matrix_row = perm->data[i];
       const int matrix_row = perm->data[i];
-      tab_text (t, 0, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[matrix_row]));
 
 
-      for (j = 0 ; j < n_factors; ++j)
+      int var_idx = pivot_category_create_leaf (
+        variables->root, pivot_value_new_variable (factor->vars[matrix_row]));
+
+      for (size_t j = 0 ; j < n_factors; ++j)
        {
          double x = gsl_matrix_get (fm, matrix_row, j);
        {
          double x = gsl_matrix_get (fm, matrix_row, j);
-
-         if ( fabs (x) < factor->blank)
+         if (fabs (x) < factor->blank)
            continue;
 
            continue;
 
-         tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL);
+          pivot_table_put2 (table, j, var_idx, pivot_value_new_number (x));
        }
     }
 
   gsl_permutation_free (perm);
 
        }
     }
 
   gsl_permutation_free (perm);
 
-  tab_submit (t);
+  pivot_table_submit (table);
 }
 
 }
 
+static void
+put_variance (struct pivot_table *table, int row, int phase_idx,
+              double lambda, double percent, double cum)
+{
+  double entries[] = { lambda, percent, cum };
+  for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
+    pivot_table_put3 (table, i, phase_idx, row,
+                      pivot_value_new_number (entries[i]));
+}
 
 static void
 
 static void
-show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
+show_explained_variance (const struct cmd_factor * factor,
+                        const struct idata *idata,
                         const gsl_vector *initial_eigenvalues,
                         const gsl_vector *extracted_eigenvalues,
                         const gsl_vector *rotated_loadings)
 {
                         const gsl_vector *initial_eigenvalues,
                         const gsl_vector *extracted_eigenvalues,
                         const gsl_vector *rotated_loadings)
 {
-  size_t i;
-  int c = 0;
-  const int heading_columns = 1;
-  const int heading_rows = 2;
-  const int nr = heading_rows + factor->n_vars;
-
-  struct tab_table *t ;
-
-  double i_total = 0.0;
-  double i_cum = 0.0;
-
-  double e_total = 0.0;
-  double e_cum = 0.0;
-
-  double r_cum = 0.0;
-
-  int nc = heading_columns;
-
-  if (factor->print & PRINT_EXTRACTION)
-    nc += 3;
-
-  if (factor->print & PRINT_INITIAL)
-    nc += 3;
-
-  if (factor->print & PRINT_ROTATION)
-    nc += 3;
-
-  /* No point having a table with only headings */
-  if ( nc <= heading_columns)
+  if (!(factor->print & (PRINT_INITIAL | PRINT_EXTRACTION | PRINT_ROTATION)))
     return;
 
     return;
 
-  t = tab_create (nc, nr);
-
-  tab_title (t, _("Total Variance Explained"));
-
-  tab_headers (t, heading_columns, 0, heading_rows, 0);
-
-  /* Outline the box */
-  tab_box (t,
-          TAL_2, TAL_2,
-          -1, -1,
-          0, 0,
-          nc - 1, nr - 1);
-
-  /* Vertical lines */
-  tab_box (t,
-          -1, -1,
-          -1, TAL_1,
-          heading_columns, 0,
-          nc - 1, nr - 1);
-
-  tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
-  tab_hline (t, TAL_1, 1, nc - 1, 1);
-
-  tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
+  struct pivot_table *table = pivot_table_create (
+    N_("Total Variance Explained"));
 
 
+  pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
+                          N_("Total"), PIVOT_RC_OTHER,
+                          /* xgettext:no-c-format */
+                          N_("% of Variance"), PIVOT_RC_PERCENT,
+                         /* xgettext:no-c-format */
+                          N_("Cumulative %"), PIVOT_RC_PERCENT);
 
 
-  if ( factor->extraction == EXTRACTION_PC)
-    tab_text (t, 0, 1, TAB_LEFT | TAT_TITLE, _("Component"));
-  else
-    tab_text (t, 0, 1, TAB_LEFT | TAT_TITLE, _("Factor"));
-
-  c = 1;
+  struct pivot_dimension *phase = pivot_dimension_create (
+    table, PIVOT_AXIS_COLUMN, N_("Phase"));
   if (factor->print & PRINT_INITIAL)
   if (factor->print & PRINT_INITIAL)
-    {
-      tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Initial Eigenvalues"));
-      c += 3;
-    }
+    pivot_category_create_leaves (phase->root, N_("Initial Eigenvalues"));
 
   if (factor->print & PRINT_EXTRACTION)
 
   if (factor->print & PRINT_EXTRACTION)
-    {
-      tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Extraction Sums of Squared Loadings"));
-      c += 3;
-    }
+    pivot_category_create_leaves (phase->root,
+                                  N_("Extraction Sums of Squared Loadings"));
 
   if (factor->print & PRINT_ROTATION)
 
   if (factor->print & PRINT_ROTATION)
-    {
-      tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Rotation Sums of Squared Loadings"));
-      c += 3;
-    }
-
-  for (i = 0; i < (nc - heading_columns) / 3 ; ++i)
-    {
-      tab_text (t, i * 3 + 1, 1, TAB_CENTER | TAT_TITLE, _("Total"));
-      /* xgettext:no-c-format */
-      tab_text (t, i * 3 + 2, 1, TAB_CENTER | TAT_TITLE, _("% of Variance"));
-      tab_text (t, i * 3 + 3, 1, TAB_CENTER | TAT_TITLE, _("Cumulative %"));
+    pivot_category_create_leaves (phase->root,
+                                  N_("Rotation Sums of Squared Loadings"));
 
 
-      tab_vline (t, TAL_2, heading_columns + i * 3, 0, nr - 1);
-    }
+  struct pivot_dimension *components = pivot_dimension_create (
+    table, PIVOT_AXIS_ROW,
+    factor->extraction == EXTRACTION_PC ? N_("Component") : N_("Factor"));
 
 
-  for (i = 0 ; i < initial_eigenvalues->size; ++i)
+  double i_total = 0.0;
+  for (size_t i = 0 ; i < initial_eigenvalues->size; ++i)
     i_total += gsl_vector_get (initial_eigenvalues, i);
 
     i_total += gsl_vector_get (initial_eigenvalues, i);
 
-  if ( factor->extraction == EXTRACTION_PAF)
-    {
-      e_total = factor->n_vars;
-    }
-  else
-    {
-      e_total = i_total;
-    }
+  double e_total = (factor->extraction == EXTRACTION_PAF
+                    ? factor->n_vars
+                    : i_total);
 
 
-  for (i = 0 ; i < factor->n_vars; ++i)
+  double i_cum = 0.0;
+  double e_cum = 0.0;
+  double r_cum = 0.0;
+  for (size_t i = 0 ; i < factor->n_vars; ++i)
     {
       const double i_lambda = gsl_vector_get (initial_eigenvalues, i);
       double i_percent = 100.0 * i_lambda / i_total ;
     {
       const double i_lambda = gsl_vector_get (initial_eigenvalues, i);
       double i_percent = 100.0 * i_lambda / i_total ;
+      i_cum += i_percent;
 
       const double e_lambda = gsl_vector_get (extracted_eigenvalues, i);
       double e_percent = 100.0 * e_lambda / e_total ;
 
       const double e_lambda = gsl_vector_get (extracted_eigenvalues, i);
       double e_percent = 100.0 * e_lambda / e_total ;
+      e_cum += e_percent;
 
 
-      const double r_lambda = gsl_vector_get (rotated_loadings, i);
-      double r_percent = 100.0 * r_lambda / e_total ;
-
-      c = 0;
+      int row = pivot_category_create_leaf (
+        components->root, pivot_value_new_integer (i + 1));
 
 
-      tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%zu"), i + 1);
-
-      i_cum += i_percent;
-      e_cum += e_percent;
-      r_cum += r_percent;
+      int phase_idx = 0;
 
       /* Initial Eigenvalues */
       if (factor->print & PRINT_INITIAL)
 
       /* Initial Eigenvalues */
       if (factor->print & PRINT_INITIAL)
-      {
-       tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL);
-       tab_double (t, c++, i + heading_rows, 0, i_percent, NULL);
-       tab_double (t, c++, i + heading_rows, 0, i_cum, NULL);
-      }
-
+        put_variance (table, row, phase_idx++, i_lambda, i_percent, i_cum);
 
 
-      if (factor->print & PRINT_EXTRACTION)
-       {
-         if (i < idata->n_extractions)
-           {
-             /* Sums of squared loadings */
-             tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL);
-             tab_double (t, c++, i + heading_rows, 0, e_percent, NULL);
-             tab_double (t, c++, i + heading_rows, 0, e_cum, NULL);
-           }
-       }
-
-      if (factor->print & PRINT_ROTATION)
-       {
-         if (i < idata->n_extractions)
-           {
-             tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL);
-             tab_double (t, c++, i + heading_rows, 0, r_percent, NULL);
-             tab_double (t, c++, i + heading_rows, 0, r_cum, NULL);
-           }
-      }
+      if (i < idata->n_extractions)
+        {
+          if (factor->print & PRINT_EXTRACTION)
+            put_variance (table, row, phase_idx++, e_lambda, e_percent, e_cum);
 
 
+          if (rotated_loadings != NULL && factor->print & PRINT_ROTATION)
+            {
+              double r_lambda = gsl_vector_get (rotated_loadings, i);
+              double r_percent = 100.0 * r_lambda / e_total ;
+              if (factor->rotation == ROT_PROMAX)
+                r_lambda = r_percent = SYSMIS;
+
+              r_cum += r_percent;
+              put_variance (table, row, phase_idx++, r_lambda, r_percent,
+                            r_cum);
+            }
+        }
     }
 
     }
 
-  tab_submit (t);
+  pivot_table_submit (table);
 }
 
 }
 
-
 static void
 static void
-show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata)
+show_factor_correlation (const struct cmd_factor * factor, const gsl_matrix *fcm)
 {
 {
-  struct tab_table *t ;
-  size_t i, j;
-  int y_pos_corr = -1;
-  int y_pos_sig = -1;
-  int suffix_rows = 0;
+  struct pivot_table *table = pivot_table_create (
+    N_("Factor Correlation Matrix"));
 
 
-  const int heading_rows = 1;
-  const int heading_columns = 2;
+  create_numeric_dimension (
+    table, PIVOT_AXIS_ROW,
+    factor->extraction == EXTRACTION_PC ? N_("Component") : N_("Factor"),
+    fcm->size2, true);
 
 
-  int nc = heading_columns ;
-  int nr = heading_rows ;
-  int n_data_sets = 0;
+  create_numeric_dimension (table, PIVOT_AXIS_COLUMN, N_("Factor 2"),
+                            fcm->size1, false);
 
 
-  if (factor->print & PRINT_CORRELATION)
-    {
-      y_pos_corr = n_data_sets;
-      n_data_sets++;
-      nc = heading_columns + factor->n_vars;
-    }
+  for (size_t i = 0 ; i < fcm->size1; ++i)
+    for (size_t j = 0 ; j < fcm->size2; ++j)
+      pivot_table_put2 (table, j, i,
+                        pivot_value_new_number (gsl_matrix_get (fcm, i, j)));
 
 
-  if (factor->print & PRINT_SIG)
-    {
-      y_pos_sig = n_data_sets;
-      n_data_sets++;
-      nc = heading_columns + factor->n_vars;
-    }
+  pivot_table_submit (table);
+}
 
 
-  nr += n_data_sets * factor->n_vars;
+static void
+add_var_dims (struct pivot_table *table, const struct cmd_factor *factor)
+{
+  for (int i = 0; i < 2; i++)
+    {
+      struct pivot_dimension *d = pivot_dimension_create (
+        table, i ? PIVOT_AXIS_ROW : PIVOT_AXIS_COLUMN,
+        N_("Variables"));
 
 
-  if (factor->print & PRINT_DETERMINANT)
-    suffix_rows = 1;
+      for (size_t j = 0; j < factor->n_vars; j++)
+        pivot_category_create_leaf (
+          d->root, pivot_value_new_variable (factor->vars[j]));
+    }
+}
 
 
-  /* If the table would contain only headings, don't bother rendering it */
-  if (nr <= heading_rows && suffix_rows == 0)
+static void
+show_aic (const struct cmd_factor *factor, const struct idata *idata)
+{
+  if ((factor->print & PRINT_AIC) == 0)
     return;
 
     return;
 
-  t = tab_create (nc, nr + suffix_rows);
-
-  tab_title (t, _("Correlation Matrix"));
-
-  tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
-
-  if (nr > heading_rows)
-    {
-      tab_headers (t, heading_columns, 0, heading_rows, 0);
+  struct pivot_table *table = pivot_table_create (N_("Anti-Image Matrices"));
 
 
-      tab_vline (t, TAL_2, 2, 0, nr - 1);
+  add_var_dims (table, factor);
 
 
-      /* Outline the box */
-      tab_box (t,
-              TAL_2, TAL_2,
-              -1, -1,
-              0, 0,
-              nc - 1, nr - 1);
+  pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Statistics"),
+                          N_("Anti-image Covariance"),
+                          N_("Anti-image Correlation"));
 
 
-      /* Vertical lines */
-      tab_box (t,
-              -1, -1,
-              -1, TAL_1,
-              heading_columns, 0,
-              nc - 1, nr - 1);
+  for (size_t i = 0; i < factor->n_vars; ++i)
+    for (size_t j = 0; j < factor->n_vars; ++j)
+      {
+        double cov = gsl_matrix_get (idata->ai_cov, i, j);
+        pivot_table_put3 (table, i, j, 0, pivot_value_new_number (cov));
 
 
+        double corr = gsl_matrix_get (idata->ai_cor, i, j);
+        pivot_table_put3 (table, i, j, 1, pivot_value_new_number (corr));
+      }
 
 
-      for (i = 0; i < factor->n_vars; ++i)
-       tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i]));
+  pivot_table_submit (table);
+}
 
 
+static void
+show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata)
+{
+  if (!(factor->print & (PRINT_CORRELATION | PRINT_SIG | PRINT_DETERMINANT)))
+    return;
 
 
-      for (i = 0 ; i < n_data_sets; ++i)
-       {
-         int y = heading_rows + i * factor->n_vars;
-         size_t v;
-         for (v = 0; v < factor->n_vars; ++v)
-           tab_text (t, 1, y + v, TAT_TITLE, var_to_string (factor->vars[v]));
+  struct pivot_table *table = pivot_table_create (N_("Correlation Matrix"));
 
 
-         tab_hline (t, TAL_1, 0, nc - 1, y);
-       }
+  if (factor->print & (PRINT_CORRELATION | PRINT_SIG))
+    {
+      add_var_dims (table, factor);
 
 
+      struct pivot_dimension *statistics = pivot_dimension_create (
+        table, PIVOT_AXIS_ROW, N_("Statistics"));
       if (factor->print & PRINT_CORRELATION)
       if (factor->print & PRINT_CORRELATION)
-       {
-         const double y = heading_rows + y_pos_corr;
-         tab_text (t, 0, y, TAT_TITLE, _("Correlations"));
-
-         for (i = 0; i < factor->n_vars; ++i)
-           {
-             for (j = 0; j < factor->n_vars; ++j)
-               tab_double (t, heading_columns + i,  y + j, 0, gsl_matrix_get (idata->corr, i, j), NULL);
-           }
-       }
-
+        pivot_category_create_leaves (statistics->root, N_("Correlation"),
+                                      PIVOT_RC_CORRELATION);
       if (factor->print & PRINT_SIG)
       if (factor->print & PRINT_SIG)
-       {
-         const double y = heading_rows + y_pos_sig * factor->n_vars;
-         tab_text (t, 0, y, TAT_TITLE, _("Sig. (1-tailed)"));
-
-         for (i = 0; i < factor->n_vars; ++i)
-           {
-             for (j = 0; j < factor->n_vars; ++j)
-               {
-                 double rho = gsl_matrix_get (idata->corr, i, j);
-                 double w = gsl_matrix_get (idata->n, i, j);
+        pivot_category_create_leaves (statistics->root, N_("Sig. (1-tailed)"),
+                                      PIVOT_RC_SIGNIFICANCE);
 
 
-                 if (i == j)
-                   continue;
+      int stat_idx = 0;
+      if (factor->print & PRINT_CORRELATION)
+        {
+          for (int i = 0; i < factor->n_vars; ++i)
+            for (int j = 0; j < factor->n_vars; ++j)
+              {
+                double corr = gsl_matrix_get (idata->mm.corr, i, j);
+                pivot_table_put3 (table, j, i, stat_idx,
+                                  pivot_value_new_number (corr));
+              }
+          stat_idx++;
+        }
 
 
-                 tab_double (t, heading_columns + i,  y + j, 0, significance_of_correlation (rho, w), NULL);
-               }
-           }
-       }
+      if (factor->print & PRINT_SIG)
+        {
+          for (int i = 0; i < factor->n_vars; ++i)
+            for (int j = 0; j < factor->n_vars; ++j)
+              if (i != j)
+                {
+                  double rho = gsl_matrix_get (idata->mm.corr, i, j);
+                  double w = gsl_matrix_get (idata->mm.n, i, j);
+                  double sig = significance_of_correlation (rho, w);
+                  pivot_table_put3 (table, j, i, stat_idx,
+                                    pivot_value_new_number (sig));
+                }
+          stat_idx++;
+        }
     }
 
   if (factor->print & PRINT_DETERMINANT)
     {
     }
 
   if (factor->print & PRINT_DETERMINANT)
     {
-      int sign = 0;
-      double det = 0.0;
+      struct pivot_value *caption = pivot_value_new_user_text_nocopy (
+        xasprintf ("%s: %.2f", _("Determinant"), idata->detR));
+      pivot_table_set_caption (table, caption);
+    }
 
 
-      const int size = idata->corr->size1;
-      gsl_permutation *p = gsl_permutation_calloc (size);
-      gsl_matrix *tmp = gsl_matrix_calloc (size, size);
-      gsl_matrix_memcpy (tmp, idata->corr);
+  pivot_table_submit (table);
+}
 
 
-      gsl_linalg_LU_decomp (tmp, p, &sign);
-      det = gsl_linalg_LU_det (tmp, sign);
-      gsl_permutation_free (p);
-      gsl_matrix_free (tmp);
+static void
+show_covariance_matrix (const struct cmd_factor *factor, const struct idata *idata)
+{
+  if (!(factor->print & PRINT_COVARIANCE))
+    return;
 
 
+  struct pivot_table *table = pivot_table_create (N_("Covariance Matrix"));
+  add_var_dims (table, factor);
 
 
-      tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant"));
-      tab_double (t, 1, nr, 0, det, NULL);
-    }
+  for (int i = 0; i < factor->n_vars; ++i)
+    for (int j = 0; j < factor->n_vars; ++j)
+      {
+        double cov = gsl_matrix_get (idata->mm.cov, i, j);
+        pivot_table_put2 (table, j, i, pivot_value_new_number (cov));
+      }
 
 
-  tab_submit (t);
+  pivot_table_submit (table);
 }
 
 
 }
 
 
-
 static void
 do_factor (const struct cmd_factor *factor, struct casereader *r)
 {
   struct ccase *c;
 static void
 do_factor (const struct cmd_factor *factor, struct casereader *r)
 {
   struct ccase *c;
-  const gsl_matrix *var_matrix;
-  const gsl_matrix *mean_matrix;
-
-  const gsl_matrix *analysis_matrix;
   struct idata *idata = idata_alloc (factor->n_vars);
 
   struct idata *idata = idata_alloc (factor->n_vars);
 
-  struct covariance *cov = covariance_1pass_create (factor->n_vars, factor->vars,
-                                             factor->wv, factor->exclude);
+  idata->cvm = covariance_1pass_create (factor->n_vars, factor->vars,
+                                       factor->wv, factor->exclude, true);
 
 
-  for ( ; (c = casereader_read (r) ); case_unref (c))
+  for (; (c = casereader_read (r)); case_unref (c))
     {
     {
-      covariance_accumulate (cov, c);
+      covariance_accumulate (idata->cvm, c);
     }
 
     }
 
-  idata->cov = covariance_calculate (cov);
+  idata->mm.cov = covariance_calculate (idata->cvm);
 
 
-  if (idata->cov == NULL)
+  if (idata->mm.cov == NULL)
     {
       msg (MW, _("The dataset contains no complete observations. No analysis will be performed."));
     {
       msg (MW, _("The dataset contains no complete observations. No analysis will be performed."));
+      covariance_destroy (idata->cvm);
       goto finish;
     }
 
       goto finish;
     }
 
-  var_matrix = covariance_moments (cov, MOMENT_VARIANCE);
-  mean_matrix = covariance_moments (cov, MOMENT_MEAN);
-  idata->n = covariance_moments (cov, MOMENT_NONE);
+  idata->mm.var_matrix = covariance_moments (idata->cvm, MOMENT_VARIANCE);
+  idata->mm.mean_matrix = covariance_moments (idata->cvm, MOMENT_MEAN);
+  idata->mm.n = covariance_moments (idata->cvm, MOMENT_NONE);
+
+  do_factor_by_matrix (factor, idata);
+
+ finish:
+  gsl_matrix_free (idata->mm.corr);
+  gsl_matrix_free (idata->mm.cov);
+
+  idata_free (idata);
+  casereader_destroy (r);
+}
 
 
-  if ( factor->method == METHOD_CORR)
+static void
+do_factor_by_matrix (const struct cmd_factor *factor, struct idata *idata)
+{
+  if (!idata->mm.cov && !(idata->mm.corr && idata->mm.var_matrix))
     {
     {
-      idata->corr = correlation_from_covariance (idata->cov, var_matrix);
-      analysis_matrix = idata->corr;
+      msg (ME, _("The dataset has no covariance matrix or a "
+                 "correlation matrix along with standard devications."));
+      return;
     }
     }
+
+  if (idata->mm.cov && !idata->mm.corr)
+    idata->mm.corr = correlation_from_covariance (idata->mm.cov, idata->mm.var_matrix);
+  if (idata->mm.corr && !idata->mm.cov)
+    idata->mm.cov = covariance_from_correlation (idata->mm.corr, idata->mm.var_matrix);
+  if (factor->method == METHOD_CORR)
+    idata->analysis_matrix = idata->mm.corr;
   else
   else
-    analysis_matrix = idata->cov;
+    idata->analysis_matrix = idata->mm.cov;
 
 
-  if ( factor->print & PRINT_UNIVARIATE)
-    {
-      const int nc = 4;
-      int i;
-      const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0;
+  gsl_matrix *r_inv;
+  r_inv  = clone_matrix (idata->mm.corr);
+  gsl_linalg_cholesky_decomp (r_inv);
+  gsl_linalg_cholesky_invert (r_inv);
 
 
+  idata->ai_cov = anti_image_cov (r_inv);
+  idata->ai_cor = anti_image_corr (r_inv, idata);
 
 
-      const int heading_columns = 1;
-      const int heading_rows = 1;
+  int i;
+  double sum_ssq_r = 0;
+  double sum_ssq_a = 0;
+  for (i = 0; i < r_inv->size1; ++i)
+    {
+      sum_ssq_r += ssq_od_n (idata->mm.corr, i);
+      sum_ssq_a += ssq_od_n (idata->ai_cor, i);
+    }
 
 
-      const int nr = heading_rows + factor->n_vars;
+  gsl_matrix_free (r_inv);
 
 
-      struct tab_table *t = tab_create (nc, nr);
-      tab_title (t, _("Descriptive Statistics"));
+  if (factor->print & PRINT_DETERMINANT
+      || factor->print & PRINT_KMO)
+    {
+      int sign = 0;
 
 
-      tab_headers (t, heading_columns, 0, heading_rows, 0);
+      const int size = idata->mm.corr->size1;
+      gsl_permutation *p = gsl_permutation_calloc (size);
+      gsl_matrix *tmp = gsl_matrix_calloc (size, size);
+      gsl_matrix_memcpy (tmp, idata->mm.corr);
 
 
-      /* Outline the box */
-      tab_box (t,
-              TAL_2, TAL_2,
-              -1, -1,
-              0, 0,
-              nc - 1, nr - 1);
+      gsl_linalg_LU_decomp (tmp, p, &sign);
+      idata->detR = gsl_linalg_LU_det (tmp, sign);
+      gsl_permutation_free (p);
+      gsl_matrix_free (tmp);
+    }
 
 
-      /* Vertical lines */
-      tab_box (t,
-              -1, -1,
-              -1, TAL_1,
-              heading_columns, 0,
-              nc - 1, nr - 1);
+  if (factor->print & PRINT_UNIVARIATE
+      && idata->mm.n && idata->mm.mean_matrix && idata->mm.var_matrix)
+    {
+      struct pivot_table *table = pivot_table_create (
+        N_("Descriptive Statistics"));
+      pivot_table_set_weight_var (table, factor->wv);
 
 
-      tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
-      tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
+      pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
+                              N_("Mean"), PIVOT_RC_OTHER,
+                              N_("Std. Deviation"), PIVOT_RC_OTHER,
+                              N_("Analysis N"), PIVOT_RC_COUNT);
 
 
-      tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
-      tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
-      tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Analysis N"));
+      struct pivot_dimension *variables = pivot_dimension_create (
+        table, PIVOT_AXIS_ROW, N_("Variables"));
 
       for (i = 0 ; i < factor->n_vars; ++i)
        {
          const struct variable *v = factor->vars[i];
 
       for (i = 0 ; i < factor->n_vars; ++i)
        {
          const struct variable *v = factor->vars[i];
-         tab_text (t, 0, i + heading_rows, TAB_LEFT | TAT_TITLE, var_to_string (v));
 
 
-         tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (mean_matrix, i, i), NULL);
-         tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (var_matrix, i, i)), NULL);
-         tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->n, i, i), wfmt);
+          int row = pivot_category_create_leaf (
+            variables->root, pivot_value_new_variable (v));
+
+          double entries[] = {
+            gsl_matrix_get (idata->mm.mean_matrix, i, i),
+            sqrt (gsl_matrix_get (idata->mm.var_matrix, i, i)),
+            gsl_matrix_get (idata->mm.n, i, i),
+          };
+          for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
+            pivot_table_put2 (table, j, row,
+                              pivot_value_new_number (entries[j]));
        }
 
        }
 
-      tab_submit (t);
+      pivot_table_submit (table);
+    }
+
+  if (factor->print & PRINT_KMO && idata->mm.n)
+    {
+      struct pivot_table *table = pivot_table_create (
+        N_("KMO and Bartlett's Test"));
+
+      struct pivot_dimension *statistics = pivot_dimension_create (
+        table, PIVOT_AXIS_ROW, N_("Statistics"),
+        N_("Kaiser-Meyer-Olkin Measure of Sampling Adequacy"), PIVOT_RC_OTHER);
+      pivot_category_create_group (
+        statistics->root, N_("Bartlett's Test of Sphericity"),
+        N_("Approx. Chi-Square"), PIVOT_RC_OTHER,
+        N_("df"), PIVOT_RC_INTEGER,
+        N_("Sig."), PIVOT_RC_SIGNIFICANCE);
+
+      /* The literature doesn't say what to do for the value of W when
+        missing values are involved.  The best thing I can think of
+        is to take the mean average. */
+      double w = 0;
+      for (i = 0; i < idata->mm.n->size1; ++i)
+       w += gsl_matrix_get (idata->mm.n, i, i);
+      w /= idata->mm.n->size1;
+
+      double xsq = ((w - 1 - (2 * factor->n_vars + 5) / 6.0)
+                    * -log (idata->detR));
+      double df = factor->n_vars * (factor->n_vars - 1) / 2;
+      double entries[] = {
+        sum_ssq_r / (sum_ssq_r + sum_ssq_a),
+        xsq,
+        df,
+        gsl_cdf_chisq_Q (xsq, df)
+      };
+      for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
+        pivot_table_put1 (table, i, pivot_value_new_number (entries[i]));
+
+      pivot_table_submit (table);
     }
 
   show_correlation_matrix (factor, idata);
     }
 
   show_correlation_matrix (factor, idata);
+  show_covariance_matrix (factor, idata);
+  if (idata->cvm)
+    covariance_destroy (idata->cvm);
 
 
-#if 1
   {
   {
+    gsl_matrix *am = matrix_dup (idata->analysis_matrix);
     gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars);
     gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars);
-    
-    gsl_eigen_symmv (matrix_dup (analysis_matrix), idata->eval, idata->evec, workspace);
+
+    gsl_eigen_symmv (am, idata->eval, idata->evec, workspace);
 
     gsl_eigen_symmv_free (workspace);
 
     gsl_eigen_symmv_free (workspace);
+    gsl_matrix_free (am);
   }
 
   gsl_eigen_symmv_sort (idata->eval, idata->evec, GSL_EIGEN_SORT_ABS_DESC);
   }
 
   gsl_eigen_symmv_sort (idata->eval, idata->evec, GSL_EIGEN_SORT_ABS_DESC);
-#endif
 
   idata->n_extractions = n_extracted_factors (factor, idata);
 
   if (idata->n_extractions == 0)
     {
 
   idata->n_extractions = n_extracted_factors (factor, idata);
 
   if (idata->n_extractions == 0)
     {
-      msg (MW, _("The FACTOR criteria result in zero factors extracted. Therefore no analysis will be performed."));
-      goto finish;
+      msg (MW, _("The %s criteria result in zero factors extracted. Therefore no analysis will be performed."), "FACTOR");
+      return;
     }
 
   if (idata->n_extractions > factor->n_vars)
     {
     }
 
   if (idata->n_extractions > factor->n_vars)
     {
-      msg (MW, _("The FACTOR criteria result in more factors than variables, which is not meaningful. No analysis will be performed."));
-      goto finish;
+      msg (MW,
+          _("The %s criteria result in more factors than variables, which is not meaningful. No analysis will be performed."),
+          "FACTOR");
+      return;
     }
     }
-    
+
   {
     gsl_matrix *rotated_factors = NULL;
   {
     gsl_matrix *rotated_factors = NULL;
+    gsl_matrix *pattern_matrix = NULL;
+    gsl_matrix *fcm = NULL;
     gsl_vector *rotated_loadings = NULL;
 
     const gsl_vector *extracted_eigenvalues = NULL;
     gsl_vector *rotated_loadings = NULL;
 
     const gsl_vector *extracted_eigenvalues = NULL;
@@ -1798,14 +2183,14 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
     struct factor_matrix_workspace *fmw = factor_matrix_workspace_alloc (idata->msr->size, idata->n_extractions);
     gsl_matrix *factor_matrix = gsl_matrix_calloc (factor->n_vars, fmw->n_factors);
 
     struct factor_matrix_workspace *fmw = factor_matrix_workspace_alloc (idata->msr->size, idata->n_extractions);
     gsl_matrix *factor_matrix = gsl_matrix_calloc (factor->n_vars, fmw->n_factors);
 
-    if ( factor->extraction == EXTRACTION_PAF)
+    if (factor->extraction == EXTRACTION_PAF)
       {
        gsl_vector *diff = gsl_vector_alloc (idata->msr->size);
       {
        gsl_vector *diff = gsl_vector_alloc (idata->msr->size);
-       struct smr_workspace *ws = ws_create (analysis_matrix);
+       struct smr_workspace *ws = ws_create (idata->analysis_matrix);
 
        for (i = 0 ; i < factor->n_vars ; ++i)
          {
 
        for (i = 0 ; i < factor->n_vars ; ++i)
          {
-           double r2 = squared_multiple_correlation (analysis_matrix, i, ws);
+           double r2 = squared_multiple_correlation (idata->analysis_matrix, i, ws);
 
            gsl_vector_set (idata->msr, i, r2);
          }
 
            gsl_vector_set (idata->msr, i, r2);
          }
@@ -1813,18 +2198,18 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
        gsl_vector_memcpy (initial_communalities, idata->msr);
 
 
        gsl_vector_memcpy (initial_communalities, idata->msr);
 
-       for (i = 0; i < factor->iterations; ++i)
+       for (i = 0; i < factor->extraction_iterations; ++i)
          {
            double min, max;
            gsl_vector_memcpy (diff, idata->msr);
 
          {
            double min, max;
            gsl_vector_memcpy (diff, idata->msr);
 
-           iterate_factor_matrix (analysis_matrix, idata->msr, factor_matrix, fmw);
-      
+           iterate_factor_matrix (idata->analysis_matrix, idata->msr, factor_matrix, fmw);
+
            gsl_vector_sub (diff, idata->msr);
 
            gsl_vector_minmax (diff, &min, &max);
            gsl_vector_sub (diff, idata->msr);
 
            gsl_vector_minmax (diff, &min, &max);
-      
-           if ( fabs (min) < factor->econverge && fabs (max) < factor->econverge)
+
+           if (fabs (min) < factor->econverge && fabs (max) < factor->econverge)
              break;
          }
        gsl_vector_free (diff);
              break;
          }
        gsl_vector_free (diff);
@@ -1841,22 +2226,28 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
        gsl_vector_memcpy (extracted_communalities, initial_communalities);
 
 
        gsl_vector_memcpy (extracted_communalities, initial_communalities);
 
-       iterate_factor_matrix (analysis_matrix, extracted_communalities, factor_matrix, fmw);
+       iterate_factor_matrix (idata->analysis_matrix, extracted_communalities, factor_matrix, fmw);
 
 
        extracted_eigenvalues = idata->eval;
       }
 
 
 
 
        extracted_eigenvalues = idata->eval;
       }
 
 
+    show_aic (factor, idata);
     show_communalities (factor, initial_communalities, extracted_communalities);
 
     show_communalities (factor, initial_communalities, extracted_communalities);
 
-
-    if ( factor->rotation != ROT_NONE)
+    if (factor->rotation != ROT_NONE)
       {
        rotated_factors = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2);
        rotated_loadings = gsl_vector_calloc (factor_matrix->size2);
       {
        rotated_factors = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2);
        rotated_loadings = gsl_vector_calloc (factor_matrix->size2);
+       if (factor->rotation == ROT_PROMAX)
+         {
+           pattern_matrix = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2);
+           fcm = gsl_matrix_calloc (factor_matrix->size2, factor_matrix->size2);
+         }
+
 
 
-       rotate (factor, factor_matrix, extracted_communalities, rotated_factors, rotated_loadings);
+       rotate (factor, factor_matrix, extracted_communalities, rotated_factors, rotated_loadings, pattern_matrix, fcm);
       }
 
     show_explained_variance (factor, idata, idata->eval, extracted_eigenvalues, rotated_loadings);
       }
 
     show_explained_variance (factor, idata, idata->eval, extracted_eigenvalues, rotated_loadings);
@@ -1866,30 +2257,41 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
     show_scree (factor, idata);
 
     show_factor_matrix (factor, idata,
     show_scree (factor, idata);
 
     show_factor_matrix (factor, idata,
-                       factor->extraction == EXTRACTION_PC ? _("Component Matrix") : _("Factor Matrix"),
+                       (factor->extraction == EXTRACTION_PC
+                         ? N_("Component Matrix") : N_("Factor Matrix")),
                        factor_matrix);
 
                        factor_matrix);
 
-    if ( factor->rotation != ROT_NONE)
+    if (factor->rotation == ROT_PROMAX)
+      {
+       show_factor_matrix (factor, idata, N_("Pattern Matrix"),
+                            pattern_matrix);
+       gsl_matrix_free (pattern_matrix);
+      }
+
+    if (factor->rotation != ROT_NONE)
       {
        show_factor_matrix (factor, idata,
       {
        show_factor_matrix (factor, idata,
-                           factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : _("Rotated Factor Matrix"),
+                           (factor->rotation == ROT_PROMAX
+                             ? N_("Structure Matrix")
+                             : factor->extraction == EXTRACTION_PC
+                             ? N_("Rotated Component Matrix")
+                            : N_("Rotated Factor Matrix")),
                            rotated_factors);
 
        gsl_matrix_free (rotated_factors);
       }
 
                            rotated_factors);
 
        gsl_matrix_free (rotated_factors);
       }
 
+    if (factor->rotation == ROT_PROMAX)
+      {
+       show_factor_correlation (factor, fcm);
+       gsl_matrix_free (fcm);
+      }
 
 
-
+    gsl_matrix_free (factor_matrix);
+    gsl_vector_free (rotated_loadings);
     gsl_vector_free (initial_communalities);
     gsl_vector_free (extracted_communalities);
   }
     gsl_vector_free (initial_communalities);
     gsl_vector_free (extracted_communalities);
   }
-
- finish:
-
-  idata_free (idata);
-
-  casereader_destroy (r);
 }
 
 
 }
 
 
-
Unnamed repository; edit this file to name it for gitweb.