X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Ffactor.c;h=a1cd333389c7276bc5a6ed7c37b776f2daf583a1;hb=445beb7fedbdee5ddaaa481e9ce77f2e2a5a2d70;hp=f3933268431c646e70ed741fd74b7856567ab011;hpb=e943de7a3e0645852c3b659889ff3ccf38dd43ec;p=pspp

diff --git a/src/language/stats/factor.c b/src/language/stats/factor.c
index f393326843..a1cd333389 100644
--- a/src/language/stats/factor.c
+++ b/src/language/stats/factor.c
@@ -1,5 +1,6 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2009, 2010, 2011, 2012 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
@@ -19,11 +20,12 @@
 #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_cdf.h>
 
+#include "data/any-reader.h"
 #include "data/casegrouper.h"
 #include "data/casereader.h"
 #include "data/casewriter.h"
@@ -35,6 +37,8 @@
 #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"
@@ -45,6 +49,7 @@
 #include "output/charts/scree.h"
 #include "output/tab.h"
 
+
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
 #define N_(msgid) msgid
@@ -87,7 +92,7 @@ enum print_opts
     PRINT_EXTRACTION  = 0x0100,
     PRINT_INITIAL     = 0x0200,
     PRINT_KMO         = 0x0400,
-    PRINT_REPR        = 0x0800, 
+    PRINT_REPR        = 0x0800,
     PRINT_FSCORE      = 0x1000
   };
 
@@ -96,6 +101,7 @@ enum rotation_type
     ROT_VARIMAX = 0,
     ROT_EQUAMAX,
     ROT_QUARTIMAX,
+    ROT_PROMAX,
     ROT_NONE
   };
 
@@ -131,14 +137,66 @@ quartimax_coefficients (double *x, double *y,
   *y = c ;
 }
 
-static const rotation_coefficients rotation_coeff[3] = {
+static const rotation_coefficients rotation_coeff[] = {
   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;
@@ -152,12 +210,14 @@ struct cmd_factor
   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;
-  int iterations;
+  int extraction_iterations;
 
   double rconverge;
 
@@ -166,13 +226,13 @@ struct cmd_factor
   bool sort;
 };
 
+
 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 */
@@ -182,6 +242,10 @@ struct idata
   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 *
@@ -204,39 +268,31 @@ idata_free (struct idata *id)
   gsl_vector_free (id->msr);
   gsl_vector_free (id->eval);
   gsl_matrix_free (id->evec);
-  if (id->cov != NULL)
-    gsl_matrix_free (id->cov);
-  if (id->corr != NULL)
-    gsl_matrix_free (CONST_CAST (gsl_matrix *, id->corr));
+  gsl_matrix_free (id->ai_cov);
+  gsl_matrix_free (id->ai_cor);
 
   free (id);
 }
 
-
-static gsl_matrix *
-anti_image (const gsl_matrix *m)
+/* 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, j;
-  gsl_matrix *a;
+  int i;
+  double ss = 0;
   assert (m->size1 == m->size2);
 
-  a = gsl_matrix_alloc (m->size1, m->size2);
-  
+  assert (j < m->size1);
+
   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);
-	}
+      if (i == j ) continue;
+      ss += pow2 (gsl_matrix_get (m, i, j));
     }
 
-  return a;
+  return ss;
 }
 
-
 /* Return the sum of all the elements excluding row N */
 static double
 ssq_od_n (const gsl_matrix *m, int n)
@@ -246,7 +302,7 @@ ssq_od_n (const gsl_matrix *m, int n)
   assert (m->size1 == m->size2);
 
   assert (n < m->size1);
-  
+
   for (i = 0; i < m->size1; ++i)
     {
       if (i == n ) continue;
@@ -260,6 +316,58 @@ ssq_od_n (const gsl_matrix *m, int n)
 }
 
 
+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
@@ -302,11 +410,11 @@ dump_vector (const gsl_vector *v)
 #endif
 
 
-static int 
+static int
 n_extracted_factors (const struct cmd_factor *factor, struct idata *idata)
 {
   int i;
-  
+
   /* If there is a cached value, then return that. */
   if ( idata->n_extractions != 0)
     return idata->n_extractions;
@@ -318,7 +426,7 @@ n_extracted_factors (const struct cmd_factor *factor, struct idata *idata)
       idata->n_extractions = factor->n_factors;
       goto finish;
     }
-  
+
   /* Use the MIN_EIGEN setting. */
   for (i = 0 ; i < idata->eval->size; ++i)
     {
@@ -353,7 +461,7 @@ struct smr_workspace
 {
   /* Copy of the subject */
   gsl_matrix *m;
-  
+
   gsl_matrix *inverse;
 
   gsl_permutation *perm;
@@ -366,7 +474,7 @@ struct smr_workspace
 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);
@@ -389,13 +497,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)
 {
-  /* 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
   */
 
@@ -407,7 +515,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);
 
-  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);
 
@@ -456,7 +564,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);
-  
+
   return ws;
 }
 
@@ -491,10 +599,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:
-   
+
    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,
@@ -527,7 +635,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);
@@ -536,7 +644,7 @@ 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);
-	  
+
 	  if ( maxindex > column_n )
 	    break;
 
@@ -554,11 +662,11 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm)
 
   gsl_permutation_free (p);
   gsl_matrix_free (mat);
-  
+
   assert ( 0 == gsl_permutation_valid (perm));
 
   /* We want the biggest value to be first */
-  gsl_permutation_reverse (perm);    
+  gsl_permutation_reverse (perm);
 }
 
 
@@ -592,7 +700,7 @@ clone_matrix (const gsl_matrix *m)
 }
 
 
-static double 
+static double
 initial_sv (const gsl_matrix *fm)
 {
   int j, k;
@@ -621,7 +729,9 @@ static void
 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;
@@ -658,7 +768,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
   /* 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)
@@ -681,7 +791,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 u = jv * jv - kv * kv;
 		  double v = 2 * jv * kv;
 		  a += u;
@@ -730,6 +840,142 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
   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*/
   for (i = 0 ; i < result->size2; ++i)
@@ -740,7 +986,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
 	{
 	  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);
@@ -761,7 +1007,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
-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;
@@ -809,12 +1055,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 void do_factor_by_matrix (const struct cmd_factor *factor, struct idata *idata);
+
+
 
 int
 cmd_factor (struct lexer *lexer, struct dataset *ds)
 {
-  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;
@@ -825,39 +1074,112 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
   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.wv = NULL;
 
   factor.rconverge = 0.0001;
 
-  factor.wv = dict_get_weight (dict);
-
   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 = create_matrix_reader_from_case_reader (dict, matrix_reader,
+						  &factor.vars, &factor.n_vars);
+    }
+  else
+    {
+      goto error;
+    }
 
   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;
+        }
+      else if (lex_match_id (lexer, "PLOT"))
 	{
           lex_match (lexer, T_EQUALS);
           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
@@ -916,6 +1238,19 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
 		{
 		  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;
@@ -926,6 +1261,7 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
 		  goto error;
 		}
 	    }
+          factor.rotation_iterations = n_iterations;
 	}
       else if (lex_match_id (lexer, "CRITERIA"))
 	{
@@ -934,59 +1270,64 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
 	    {
 	      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);
-		      lex_force_match (lexer, T_RPAREN);
+		      if (! lex_force_match (lexer, T_RPAREN))
+			goto error;
 		    }
 		}
 	      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);
-		      lex_force_match (lexer, T_RPAREN);
+		      if (! lex_force_match (lexer, T_RPAREN))
+			goto error;
 		    }
 		}
 	      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);
-		      lex_force_match (lexer, T_RPAREN);
+		      if (! lex_force_match (lexer, T_RPAREN))
+			goto error;
 		    }
 		}
 	      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"))
 		{
-		  if ( lex_force_match (lexer, T_LPAREN))
+		  if ( lex_force_match (lexer, T_LPAREN)
+                       && lex_force_int (lexer))
 		    {
-		      lex_force_int (lexer);
-		      factor.iterations = lex_integer (lexer);
+		      n_iterations = 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, "DEFAULT"))
 		{
 		  factor.n_factors = 0;
 		  factor.min_eigen = 1;
-		  factor.iterations = 25;
+		  n_iterations = 25;
 		}
 	      else
 		{
@@ -1022,6 +1363,7 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
 		  goto error;
 		}
 	    }
+          factor.extraction_iterations = n_iterations;
 	}
       else if (lex_match_id (lexer, "FORMAT"))
 	{
@@ -1034,12 +1376,13 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
 		}
 	      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);
-		      lex_force_match (lexer, T_RPAREN);
+		      if (! lex_force_match (lexer, T_RPAREN))
+			goto error;
 		    }
 		}
 	      else if (lex_match_id (lexer, "DEFAULT"))
@@ -1072,10 +1415,11 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
 	      else if (lex_match_id (lexer, "INV"))
 		{
 		}
+#endif
 	      else if (lex_match_id (lexer, "AIC"))
 		{
+		  factor.print |= PRINT_AIC;
 		}
-#endif
 	      else if (lex_match_id (lexer, "SIG"))
 		{
 		  factor.print |= PRINT_SIG;
@@ -1084,11 +1428,10 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
 		{
 		  factor.print |= PRINT_CORRELATION;
 		}
-#if FACTOR_FULLY_IMPLEMENTED
 	      else if (lex_match_id (lexer, "COVARIANCE"))
 		{
+		  factor.print |= PRINT_COVARIANCE;
 		}
-#endif
 	      else if (lex_match_id (lexer, "ROTATION"))
 		{
 		  factor.print |= PRINT_ROTATION;
@@ -1172,13 +1515,43 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
   if ( factor.rotation == ROT_NONE )
     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 (next_matrix_from_reader (&id->mm, mr,
+				      factor.vars, factor.n_vars))
+	{
+	  do_factor_by_matrix (&factor, id);
 
+	  gsl_matrix_free (id->mm.corr);
+	  id->mm.corr = NULL;
+	  gsl_matrix_free (id->mm.cov);
+	  id->mm.cov = NULL;
+	}
+
+      idata_free (id);
+    }
+  else
+    if ( ! run_factor (ds, &factor))
+      goto error;
+
+
+  destroy_matrix_reader (mr);
   free (factor.vars);
   return CMD_SUCCESS;
 
  error:
+  destroy_matrix_reader (mr);
   free (factor.vars);
   return CMD_FAILURE;
 }
@@ -1238,14 +1611,14 @@ 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
-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
-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 ;
@@ -1319,10 +1692,10 @@ show_communalities (const struct cmd_factor * factor,
       tab_text (t, c++, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[i]));
 
       if (factor->print & PRINT_INITIAL)
-	tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL);
+	tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL, RC_OTHER);
 
       if (factor->print & PRINT_EXTRACTION)
-	tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL);
+	tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL, RC_OTHER);
     }
 
   tab_submit (t);
@@ -1330,9 +1703,10 @@ show_communalities (const struct cmd_factor * factor,
 
 
 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;
@@ -1343,10 +1717,10 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const
 
   struct tab_table *t = tab_create (nc, nr);
 
-  /* 
+  /*
   if ( factor->extraction == EXTRACTION_PC )
     tab_title (t, _("Component Matrix"));
-  else 
+  else
     tab_title (t, _("Factor Matrix"));
   */
 
@@ -1411,7 +1785,7 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const
 	  if ( fabs (x) < factor->blank)
 	    continue;
 
-	  tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL);
+	  tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL, RC_OTHER);
 	}
     }
 
@@ -1422,7 +1796,8 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const
 
 
 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)
@@ -1452,7 +1827,9 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
     nc += 3;
 
   if (factor->print & PRINT_ROTATION)
-    nc += 3;
+    {
+      nc += factor->rotation == ROT_PROMAX ? 1 : 3;
+    }
 
   /* No point having a table with only headings */
   if ( nc <= heading_columns)
@@ -1504,18 +1881,23 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
 
   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;
+      const int width = factor->rotation == ROT_PROMAX ? 0 : 2;
+      tab_joint_text (t, c, 0, c + width, 0, TAB_CENTER | TAT_TITLE, _("Rotation Sums of Squared Loadings"));
+      c += width + 1;
     }
 
-  for (i = 0; i < (nc - heading_columns) / 3 ; ++i)
+  for (i = 0; i < (nc - heading_columns + 2) / 3 ; ++i)
     {
       tab_text (t, i * 3 + 1, 1, TAB_CENTER | TAT_TITLE, _("Total"));
+
+      tab_vline (t, TAL_2, heading_columns + i * 3, 0, nr - 1);
+
+      if (i == 2 && factor->rotation == ROT_PROMAX)
+	continue;
+
       /* 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 %"));
-
-      tab_vline (t, TAL_2, heading_columns + i * 3, 0, nr - 1);
     }
 
   for (i = 0 ; i < initial_eigenvalues->size; ++i)
@@ -1548,9 +1930,9 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
       /* 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);
+	tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL, RC_OTHER);
+	tab_double (t, c++, i + heading_rows, 0, i_percent, NULL, RC_OTHER);
+	tab_double (t, c++, i + heading_rows, 0, i_cum, NULL, RC_OTHER);
       }
 
 
@@ -1559,9 +1941,9 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
 	  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);
+	      tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL, RC_OTHER);
+	      tab_double (t, c++, i + heading_rows, 0, e_percent, NULL, RC_OTHER);
+	      tab_double (t, c++, i + heading_rows, 0, e_cum, NULL, RC_OTHER);
 	    }
 	}
 
@@ -1575,9 +1957,12 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
               if (i < idata->n_extractions)
                 {
                   r_cum += r_percent;
-                  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);
+                  tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL, RC_OTHER);
+		  if (factor->rotation != ROT_PROMAX)
+		    {
+		      tab_double (t, c++, i + heading_rows, 0, r_percent, NULL, RC_OTHER);
+		      tab_double (t, c++, i + heading_rows, 0, r_cum, NULL, RC_OTHER);
+		    }
                 }
             }
         }
@@ -1587,6 +1972,142 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
 }
 
 
+static void
+show_factor_correlation (const struct cmd_factor * factor, const gsl_matrix *fcm)
+{
+  size_t i, j;
+  const int heading_columns = 1;
+  const int heading_rows = 1;
+  const int nr = heading_rows + fcm->size2;
+  const int nc = heading_columns + fcm->size1;
+  struct tab_table *t = tab_create (nc, nr);
+
+  tab_title (t, _("Factor Correlation Matrix"));
+
+  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);
+
+
+  if ( factor->extraction == EXTRACTION_PC)
+    tab_text (t, 0, 0, TAB_LEFT | TAT_TITLE, _("Component"));
+  else
+    tab_text (t, 0, 0, TAB_LEFT | TAT_TITLE, _("Factor"));
+
+  for (i = 0 ; i < fcm->size1; ++i)
+    {
+      tab_text_format (t, heading_columns + i, 0, TAB_CENTER | TAT_TITLE, _("%zu"), i + 1);
+    }
+
+  for (i = 0 ; i < fcm->size2; ++i)
+    {
+      tab_text_format (t, 0, heading_rows + i, TAB_CENTER | TAT_TITLE, _("%zu"), i + 1);
+    }
+
+
+  for (i = 0 ; i < fcm->size1; ++i)
+    {
+      for (j = 0 ; j < fcm->size2; ++j)
+	tab_double (t, heading_columns + j,  heading_rows + i, 0,
+		    gsl_matrix_get (fcm, i, j), NULL, RC_OTHER);
+    }
+
+  tab_submit (t);
+}
+
+static void
+show_aic (const struct cmd_factor *factor, const struct idata *idata)
+{
+  struct tab_table *t ;
+  size_t i;
+
+  const int heading_rows = 1;
+  const int heading_columns = 2;
+
+  const int nc = heading_columns + factor->n_vars;
+  const int nr = heading_rows + 2 * factor->n_vars;
+
+  if ((factor->print & PRINT_AIC) == 0)
+    return;
+
+  t = tab_create (nc, nr);
+
+  tab_title (t, _("Anti-Image Matrices"));
+
+  tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
+
+  tab_headers (t, heading_columns, 0, heading_rows, 0);
+
+  tab_vline (t, TAL_2, 2, 0, nr - 1);
+
+  /* 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);
+
+
+  for (i = 0; i < factor->n_vars; ++i)
+    tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i]));
+
+  tab_text (t, 0, heading_rows, TAT_TITLE, _("Anti-image Covariance"));
+  tab_hline (t, TAL_1, 0, nc - 1, heading_rows + factor->n_vars);
+  tab_text (t, 0, heading_rows + factor->n_vars, TAT_TITLE, _("Anti-image Correlation"));
+
+  for (i = 0; i < factor->n_vars; ++i)
+    {
+      tab_text (t, 1, i + heading_rows, TAT_TITLE,
+		var_to_string (factor->vars[i]));
+
+      tab_text (t, 1, factor->n_vars + i + heading_rows, TAT_TITLE,
+		var_to_string (factor->vars[i]));
+    }
+
+  for (i = 0; i < factor->n_vars; ++i)
+    {
+      int j;
+      for (j = 0; j < factor->n_vars; ++j)
+	{
+	  tab_double (t, heading_columns + i, heading_rows + j, 0,
+		      gsl_matrix_get (idata->ai_cov, i, j), NULL, RC_OTHER);
+	}
+
+
+      for (j = 0; j < factor->n_vars; ++j)
+	{
+	  tab_double (t, heading_columns + i, factor->n_vars + heading_rows + j, 0,
+		      gsl_matrix_get (idata->ai_cor, i, j), NULL, RC_OTHER);
+	}
+    }
+
+  tab_submit (t);
+}
+
 static void
 show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata)
 {
@@ -1675,7 +2196,7 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id
 	  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);
+		tab_double (t, heading_columns + j,  y + i, 0, gsl_matrix_get (idata->mm.corr, i, j), NULL, RC_OTHER);
 	    }
 	}
 
@@ -1688,13 +2209,13 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id
 	    {
 	      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);
+		  double rho = gsl_matrix_get (idata->mm.corr, i, j);
+		  double w = gsl_matrix_get (idata->mm.n, i, j);
 
 		  if (i == j)
 		    continue;
 
-		  tab_double (t, heading_columns + i,  y + j, 0, significance_of_correlation (rho, w), NULL);
+		  tab_double (t, heading_columns + j,  y + i, 0, significance_of_correlation (rho, w), NULL, RC_PVALUE);
 		}
 	    }
 	}
@@ -1704,65 +2225,180 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id
     {
       tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant"));
 
-      tab_double (t, 1, nr, 0, idata->detR, NULL);
+      tab_double (t, 1, nr, 0, idata->detR, NULL, RC_OTHER);
     }
 
   tab_submit (t);
 }
 
+static void
+show_covariance_matrix (const struct cmd_factor *factor, const struct idata *idata)
+{
+  struct tab_table *t ;
+  size_t i, j;
+  int y_pos_corr = -1;
+  int suffix_rows = 0;
+
+  const int heading_rows = 1;
+  const int heading_columns = 1;
+
+  int nc = heading_columns ;
+  int nr = heading_rows ;
+  int n_data_sets = 0;
+
+  if (factor->print & PRINT_COVARIANCE)
+    {
+      y_pos_corr = n_data_sets;
+      n_data_sets++;
+      nc = heading_columns + factor->n_vars;
+    }
+
+  nr += n_data_sets * factor->n_vars;
+
+  /* If the table would contain only headings, don't bother rendering it */
+  if (nr <= heading_rows && suffix_rows == 0)
+    return;
+
+  t = tab_create (nc, nr + suffix_rows);
+
+  tab_title (t, _("Covariance Matrix"));
+
+  tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
+
+  if (nr > heading_rows)
+    {
+      tab_headers (t, heading_columns, 0, heading_rows, 0);
+
+      tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
+
+      /* 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);
+
+
+      for (i = 0; i < factor->n_vars; ++i)
+	tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i]));
+
+
+      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, heading_columns -1, y + v, TAT_TITLE, var_to_string (factor->vars[v]));
+
+	  tab_hline (t, TAL_1, 0, nc - 1, y);
+	}
+
+      if (factor->print & PRINT_COVARIANCE)
+	{
+	  const double y = heading_rows + y_pos_corr;
+
+	  for (i = 0; i < factor->n_vars; ++i)
+	    {
+	      for (j = 0; j < factor->n_vars; ++j)
+		tab_double (t, heading_columns + j,  y + i, 0, gsl_matrix_get (idata->mm.cov, i, j), NULL, RC_OTHER);
+	    }
+	}
+    }
+
+  tab_submit (t);
+}
 
 
 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 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))
     {
-      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."));
-      covariance_destroy (cov);
+      covariance_destroy (idata->cvm);
       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);
 
-  if ( factor->method == METHOD_CORR)
+  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);
+}
+
+static void
+do_factor_by_matrix (const struct cmd_factor *factor, struct idata *idata)
+{
+  if (!idata->mm.cov && !idata->mm.corr)
     {
-      idata->corr = correlation_from_covariance (idata->cov, var_matrix);
-      
-      analysis_matrix = idata->corr;
+      msg (ME, _("The dataset has no complete covariance or correlation matrix."));
+      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
-    analysis_matrix = idata->cov;
+    idata->analysis_matrix = idata->mm.cov;
+
+  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);
 
+  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 (r_inv, i);
+      sum_ssq_a += ssq_od_n (idata->ai_cov, i);
+    }
+
+  gsl_matrix_free (r_inv);
 
   if (factor->print & PRINT_DETERMINANT
       || factor->print & PRINT_KMO)
     {
       int sign = 0;
 
-      const int size = idata->corr->size1;
+      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->corr);
+      gsl_matrix_memcpy (tmp, idata->mm.corr);
 
       gsl_linalg_LU_decomp (tmp, p, &sign);
       idata->detR = gsl_linalg_LU_det (tmp, sign);
@@ -1782,6 +2418,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
       const int nr = heading_rows + factor->n_vars;
 
       struct tab_table *t = tab_create (nc, nr);
+      tab_set_format (t, RC_WEIGHT, wfmt);
       tab_title (t, _("Descriptive Statistics"));
 
       tab_headers (t, heading_columns, 0, heading_rows, 0);
@@ -1812,9 +2449,9 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 	  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);
+	  tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (idata->mm.mean_matrix, i, i), NULL, RC_OTHER);
+	  tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (idata->mm.var_matrix, i, i)), NULL, RC_OTHER);
+	  tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->mm.n, i, i), NULL, RC_WEIGHT);
 	}
 
       tab_submit (t);
@@ -1823,10 +2460,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
   if (factor->print & PRINT_KMO)
     {
       int i;
-      double sum_ssq_r = 0;
-      double sum_ssq_a = 0;
-
-      double df = factor->n_vars * ( factor->n_vars - 1) / 2;
+      double df = factor->n_vars * (factor->n_vars - 1) / 2;
 
       double w = 0;
 
@@ -1839,25 +2473,11 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
       const int nr = heading_rows + 4;
       const int nc = heading_columns + 1;
 
-      gsl_matrix *a, *x;
+
 
       struct tab_table *t = tab_create (nc, nr);
       tab_title (t, _("KMO and Bartlett's Test"));
 
-      x  = clone_matrix (idata->corr);
-      gsl_linalg_cholesky_decomp (x);
-      gsl_linalg_cholesky_invert (x);
-
-      a = anti_image (x);
-
-      for (i = 0; i < x->size1; ++i)
-	{
-	  sum_ssq_r += ssq_od_n (x, i);
-	  sum_ssq_a += ssq_od_n (a, i);
-	}
-
-      gsl_matrix_free (a);
-      gsl_matrix_free (x);
 
       tab_headers (t, heading_columns, 0, heading_rows, 0);
 
@@ -1872,7 +2492,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
       tab_text (t, 0, 0, TAT_TITLE | TAB_LEFT, _("Kaiser-Meyer-Olkin Measure of Sampling Adequacy"));
 
-      tab_double (t, 2, 0, 0, sum_ssq_r /  (sum_ssq_r + sum_ssq_a), NULL);
+      tab_double (t, 2, 0, 0, sum_ssq_r /  (sum_ssq_r + sum_ssq_a), NULL, RC_OTHER);
 
       tab_text (t, 0, 1, TAT_TITLE | TAB_LEFT, _("Bartlett's Test of Sphericity"));
 
@@ -1881,32 +2501,34 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
       tab_text (t, 1, 3, TAT_TITLE, _("Sig."));
 
 
-      /* The literature doesn't say what to do for the value of W when 
+      /* 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. */
       w = 0;
-      for (i = 0; i < idata->n->size1; ++i)
-	w += gsl_matrix_get (idata->n, i, i);
-      w /= idata->n->size1;
+      for (i = 0; i < idata->mm.n->size1; ++i)
+	w += gsl_matrix_get (idata->mm.n, i, i);
+      w /= idata->mm.n->size1;
 
       xsq = w - 1 - (2 * factor->n_vars + 5) / 6.0;
       xsq *= -log (idata->detR);
 
-      tab_double (t, 2, 1, 0, xsq, NULL);
-      tab_double (t, 2, 2, 0, df, &F_8_0);
-      tab_double (t, 2, 3, 0, gsl_cdf_chisq_Q (xsq, df), NULL);
-      
+      tab_double (t, 2, 1, 0, xsq, NULL, RC_OTHER);
+      tab_double (t, 2, 2, 0, df, NULL, RC_INTEGER);
+      tab_double (t, 2, 3, 0, gsl_cdf_chisq_Q (xsq, df), NULL, RC_PVALUE);
+
 
       tab_submit (t);
     }
 
   show_correlation_matrix (factor, idata);
-  covariance_destroy (cov);
+  show_covariance_matrix (factor, idata);
+  if (idata->cvm)
+    covariance_destroy (idata->cvm);
 
   {
-    gsl_matrix *am = matrix_dup (analysis_matrix);
+    gsl_matrix *am = matrix_dup (idata->analysis_matrix);
     gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars);
-    
+
     gsl_eigen_symmv (am, idata->eval, idata->evec, workspace);
 
     gsl_eigen_symmv_free (workspace);
@@ -1919,18 +2541,22 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
   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)
     {
-      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 *pattern_matrix = NULL;
+    gsl_matrix *fcm = NULL;
     gsl_vector *rotated_loadings = NULL;
 
     const gsl_vector *extracted_eigenvalues = NULL;
@@ -1943,11 +2569,11 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
     if ( factor->extraction == EXTRACTION_PAF)
       {
 	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)
 	  {
-	    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);
 	  }
@@ -1955,17 +2581,17 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
 	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);
 
-	    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);
-      
+
 	    if ( fabs (min) < factor->econverge && fabs (max) < factor->econverge)
 	      break;
 	  }
@@ -1983,22 +2609,28 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
 	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;
       }
 
 
+    show_aic (factor, idata);
     show_communalities (factor, initial_communalities, extracted_communalities);
 
-
     if ( factor->rotation != ROT_NONE)
       {
 	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);
@@ -2011,29 +2643,34 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 			factor->extraction == EXTRACTION_PC ? _("Component Matrix") : _("Factor Matrix"),
 			factor_matrix);
 
+    if ( factor->rotation == ROT_PROMAX)
+      {
+	show_factor_matrix (factor, idata, _("Pattern Matrix"),  pattern_matrix);
+	gsl_matrix_free (pattern_matrix);
+      }
+
     if ( factor->rotation != ROT_NONE)
       {
 	show_factor_matrix (factor, idata,
-			    factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : _("Rotated Factor Matrix"),
+			    (factor->rotation == ROT_PROMAX) ? _("Structure Matrix") :
+			    (factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") :
+			     _("Rotated Factor Matrix")),
 			    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);
   }
-
- finish:
-
-  idata_free (idata);
-
-  casereader_destroy (r);
 }
 
 
-