diff --git a/compute/zgesvd.c b/compute/zgesvd.c
index 65954ded796438f509a8f66957588cf130ee4312..7d12a10ba10f1f38c23f00501a5bb424bfaf6c8b 100644
--- a/compute/zgesvd.c
+++ b/compute/zgesvd.c
@@ -44,7 +44,7 @@
  *  are returned in descending order.  The first min(m,n) columns of
  *  U and V are the left and right singular vectors of A.
  *
- *  Note that the routine returns V**T, not V.
+ *  Note that the routine returns V^T, not V.
  *******************************************************************************
  *
  * @param[in] jobu
diff --git a/compute/zgetrs_incpiv.c b/compute/zgetrs_incpiv.c
index 092066b6ea264389674837d2fc18f91f696547bc..225cb125f7c748a49f0816d1d21ca1985df6e361 100644
--- a/compute/zgetrs_incpiv.c
+++ b/compute/zgetrs_incpiv.c
@@ -37,7 +37,7 @@
  * @param[in] trans
  *          Intended to specify the the form of the system of equations:
  *          = ChamNoTrans:   A * X = B     (No transpose)
- *          = ChamTrans:     A**T * X = B  (Transpose)
+ *          = ChamTrans:     A^T * X = B  (Transpose)
  *          = ChamConjTrans: A^H * X = B  (Conjugate transpose)
  *          Currently only ChamNoTrans is supported.
  *
diff --git a/compute/zgetrs_nopiv.c b/compute/zgetrs_nopiv.c
index 9392d31a1b723f79b7e596d50e8506a522b85697..3a3dfe36041d3a8ba4c6f300a4d40a3072578a46 100644
--- a/compute/zgetrs_nopiv.c
+++ b/compute/zgetrs_nopiv.c
@@ -38,7 +38,7 @@
  * @param[in] trans
  *          Intended to specify the the form of the system of equations:
  *          = ChamNoTrans:   A * X = B     (No transpose)
- *          = ChamTrans:     A**T * X = B  (Transpose)
+ *          = ChamTrans:     A^T * X = B  (Transpose)
  *          = ChamConjTrans: A^H * X = B  (Conjugate transpose)
  *          Currently only ChamNoTrans is supported.
  *
diff --git a/compute/zsysv.c b/compute/zsysv.c
index 256e27de9f7a61b92ee0243b50d792aa7eabb8c0..ebee2c6f13de6b154a371c1f04bb9fc484a1c019 100644
--- a/compute/zsysv.c
+++ b/compute/zsysv.c
@@ -62,7 +62,7 @@
  *          triangular part of the matrix A, and the strictly upper triangular part of A is not
  *          referenced.
  *          On exit, if return value = 0, the factor U or L from the Cholesky factorization
- *          A = U**T*U or A = L*L**T.
+ *          A = U^T*U or A = L*L^T.
  *
  * @param[in] LDA
  *          The leading dimension of the array A. LDA >= max(1,N).
@@ -194,7 +194,7 @@ int CHAMELEON_zsysv( cham_uplo_t uplo, int N, int NRHS,
  *          triangular part of the matrix A, and the strictly upper triangular part of A is not
  *          referenced.
  *          On exit, if return value = 0, the factor U or L from the Cholesky factorization
- *          A = U**T*U or A = L*L**T.
+ *          A = U^T*U or A = L*L^T.
  *
  * @param[in,out] B
  *          On entry, the N-by-NRHS right hand side matrix B.
diff --git a/compute/zsytrf.c b/compute/zsytrf.c
index 3f13655cec08418d4fb60387ac5eeef3129065d5..b603ddde5224e523317fc9580eb9eb87b1704080 100644
--- a/compute/zsytrf.c
+++ b/compute/zsytrf.c
@@ -164,7 +164,7 @@ int CHAMELEON_zsytrf( cham_uplo_t uplo, int N,
  *          triangular part of the matrix A, and the strictly upper triangular part of A is not
  *          referenced.
  *          On exit, if return value = 0, the factor U or L from the Cholesky factorization
- *          A = U**T*U or A = L*L**T.
+ *          A = U^T*U or A = L*L^T.
  *
  *******************************************************************************
  *
diff --git a/compute/zsytrs.c b/compute/zsytrs.c
index 7cf60f28e746dcaf95de3a531840ad7e60f3568c..84a2c778fc1bcf36bfe843a6194ec6cab27ad373 100644
--- a/compute/zsytrs.c
+++ b/compute/zsytrs.c
@@ -48,7 +48,7 @@
  *          The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.
  *
  * @param[in] A
- *          The triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T,
+ *          The triangular factor U or L from the Cholesky factorization A = U^T*U or A = L*L^T,
  *          computed by CHAMELEON_zsytrf.
  *
  * @param[in] LDA
@@ -172,7 +172,7 @@ int CHAMELEON_zsytrs( cham_uplo_t uplo, int N, int NRHS,
  *          = ChamLower: Lower triangle of A is stored.
  *
  * @param[in] A
- *          The triangular factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T,
+ *          The triangular factor U or L from the Cholesky factorization A = U^T*U or A = L*L^T,
  *          computed by CHAMELEON_zsytrf.
  *
  * @param[in,out] B
diff --git a/coreblas/compute/core_zherfb.c b/coreblas/compute/core_zherfb.c
index d71c18028e8ed0c24c2e3b57e50a9c8a69b15955..a54de4a5aaec39def2aa98b11ae250bc8dff96f4 100644
--- a/coreblas/compute/core_zherfb.c
+++ b/coreblas/compute/core_zherfb.c
@@ -27,7 +27,7 @@
  *
  *  CORE_zherfb overwrites the symmetric complex N-by-N tile C with
  *
- *    Q**T*C*Q
+ *    Q^T*C*Q
  *
  *  where Q is a complex unitary matrix defined as the product of k
  *  elementary reflectors
@@ -72,7 +72,7 @@
  *
  * @param[in,out] C
  *         On entry, the symmetric N-by-N tile C.
- *         On exit, C is overwritten by Q**T*C*Q.
+ *         On exit, C is overwritten by Q^T*C*Q.
  *
  * @param[in] ldc
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/coreblas/compute/core_zpamm.c b/coreblas/compute/core_zpamm.c
index 438737436612e02277f98926adee1cb5fd5ed216..01a25ea5556a8e959d8a9e16befe070fa8ef9359 100644
--- a/coreblas/compute/core_zpamm.c
+++ b/coreblas/compute/core_zpamm.c
@@ -52,7 +52,7 @@ static inline int CORE_zpamm_w(cham_side_t side, cham_trans_t trans, cham_uplo_t
  *
  *  where  op( V ) is one of
  *
- *     op( V ) = V   or   op( V ) = V**T   or   op( V ) = V^H,
+ *     op( V ) = V   or   op( V ) = V^T   or   op( V ) = V^H,
  *
  *  A1, A2 and W are general matrices, and V is:
  *
diff --git a/coreblas/compute/core_zpemv.c b/coreblas/compute/core_zpemv.c
index 846c2695078d88a630e97db4ad05827970a1c7e5..0144020a748f0ab3a2368fa819411a80aa6441b7 100644
--- a/coreblas/compute/core_zpemv.c
+++ b/coreblas/compute/core_zpemv.c
@@ -36,7 +36,7 @@
  *
  *  where  op( A ) is one of
  *
- *     op( A ) = A   or   op( A ) = A**T   or   op( A ) = A^H,
+ *     op( A ) = A   or   op( A ) = A^T   or   op( A ) = A^H,
  *
  *  alpha and beta are scalars, x and y are vectors and A is a
  *  pentagonal matrix (see further details).
@@ -52,7 +52,7 @@
  * @param[in] trans
  *
  *         @arg ChamNoTrans   :  y := alpha*A*x    + beta*y.
- *         @arg ChamTrans     :  y := alpha*A**T*x + beta*y.
+ *         @arg ChamTrans     :  y := alpha*A^T*x + beta*y.
  *         @arg ChamConjTrans :  y := alpha*A^H*x + beta*y.
  *
  * @param[in] M
diff --git a/coreblas/compute/core_zunmlq.c b/coreblas/compute/core_zunmlq.c
index b1b98c4853f9ec8bb9bf2b480f986876369c0cf2..6a310b8598349a5c6a1adf9d2c086ef5a13548eb 100644
--- a/coreblas/compute/core_zunmlq.c
+++ b/coreblas/compute/core_zunmlq.c
@@ -90,7 +90,7 @@
  *
  * @param[in,out] C
  *         On entry, the M-by-N tile C.
- *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+ *         On exit, C is overwritten by Q*C or Q^T*C or C*Q^T or C*Q.
  *
  * @param[in] LDC
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/coreblas/compute/core_zunmqr.c b/coreblas/compute/core_zunmqr.c
index 068a5f40c031330c3fc8ff8120784b51b9ee5841..712da7e6c11a0c6accc01ca2418446eebb379637 100644
--- a/coreblas/compute/core_zunmqr.c
+++ b/coreblas/compute/core_zunmqr.c
@@ -91,7 +91,7 @@
  *
  * @param[in,out] C
  *         On entry, the M-by-N tile C.
- *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+ *         On exit, C is overwritten by Q*C or Q^T*C or C*Q^T or C*Q.
  *
  * @param[in] LDC
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/runtime/openmp/codelets/codelet_zunmlq.c b/runtime/openmp/codelets/codelet_zunmlq.c
index 6a4f7bdac9fe8456e73f0551f305560b84d4eab6..92d6e71f8f34c171100880d747fb4f605e50baa3 100644
--- a/runtime/openmp/codelets/codelet_zunmlq.c
+++ b/runtime/openmp/codelets/codelet_zunmlq.c
@@ -90,7 +90,7 @@
  *
  * @param[in,out] C
  *         On entry, the M-by-N tile C.
- *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+ *         On exit, C is overwritten by Q*C or Q^T*C or C*Q^T or C*Q.
  *
  * @param[in] LDC
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/runtime/openmp/codelets/codelet_zunmqr.c b/runtime/openmp/codelets/codelet_zunmqr.c
index 2e796545ff105355f9b854b9fab426f4697dae41..66aa62b5dbc6e8ce2d9a74b91ce1fd5ae1d088f8 100644
--- a/runtime/openmp/codelets/codelet_zunmqr.c
+++ b/runtime/openmp/codelets/codelet_zunmqr.c
@@ -90,7 +90,7 @@
  *
  * @param[in,out] C
  *         On entry, the M-by-N tile C.
- *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+ *         On exit, C is overwritten by Q*C or Q^T*C or C*Q^T or C*Q.
  *
  * @param[in] LDC
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/runtime/quark/codelets/codelet_zunmlq.c b/runtime/quark/codelets/codelet_zunmlq.c
index bed3fa03fe52cf8b1b2ff0b981e4d845ee8d2e78..5b8687571af808182585b223ed9636099fccc296 100644
--- a/runtime/quark/codelets/codelet_zunmlq.c
+++ b/runtime/quark/codelets/codelet_zunmlq.c
@@ -114,7 +114,7 @@ void CORE_zunmlq_quark(Quark *quark)
  *
  * @param[in,out] C
  *         On entry, the M-by-N tile C.
- *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+ *         On exit, C is overwritten by Q*C or Q^T*C or C*Q^T or C*Q.
  *
  * @param[in] LDC
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/runtime/quark/codelets/codelet_zunmqr.c b/runtime/quark/codelets/codelet_zunmqr.c
index 34ca5fdce9751c400610b2ac08ee808b4b23b882..f03746016c3012467e4a1370353084f8cb21f282 100644
--- a/runtime/quark/codelets/codelet_zunmqr.c
+++ b/runtime/quark/codelets/codelet_zunmqr.c
@@ -114,7 +114,7 @@ void CORE_zunmqr_quark(Quark *quark)
  *
  * @param[in,out] C
  *         On entry, the M-by-N tile C.
- *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+ *         On exit, C is overwritten by Q*C or Q^T*C or C*Q^T or C*Q.
  *
  * @param[in] LDC
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/runtime/starpu/codelets/codelet_zunmlq.c b/runtime/starpu/codelets/codelet_zunmlq.c
index 28d4c509549fc50f0c405cb797691c50142ddd26..be36f957d6451521ca38415dfd04eaba3e332e0b 100644
--- a/runtime/starpu/codelets/codelet_zunmlq.c
+++ b/runtime/starpu/codelets/codelet_zunmlq.c
@@ -159,7 +159,7 @@ CODELETS(zunmlq, 4, cl_zunmlq_cpu_func, cl_zunmlq_cuda_func, STARPU_CUDA_ASYNC)
  *
  * @param[in,out] C
  *         On entry, the M-by-N tile C.
- *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+ *         On exit, C is overwritten by Q*C or Q^T*C or C*Q^T or C*Q.
  *
  * @param[in] LDC
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/runtime/starpu/codelets/codelet_zunmqr.c b/runtime/starpu/codelets/codelet_zunmqr.c
index 6681a39f74c9a109f5b0e6d8c2ed3469d1066933..8ff98bc7925a5d360d75f2f573fc5dd0674f1abb 100644
--- a/runtime/starpu/codelets/codelet_zunmqr.c
+++ b/runtime/starpu/codelets/codelet_zunmqr.c
@@ -159,7 +159,7 @@ CODELETS(zunmqr, 4, cl_zunmqr_cpu_func, cl_zunmqr_cuda_func, STARPU_CUDA_ASYNC)
  *
  * @param[in,out] C
  *         On entry, the M-by-N tile C.
- *         On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
+ *         On exit, C is overwritten by Q*C or Q^T*C or C*Q^T or C*Q.
  *
  * @param[in] LDC
  *         The leading dimension of the array C. LDC >= max(1,M).
diff --git a/testing/lin/clagsy.f b/testing/lin/clagsy.f
index c5fea1b7d09103272096213443dc251191f2c3fb..0522d060097c4d5c6790857a684f8a8067b0dfc1 100644
--- a/testing/lin/clagsy.f
+++ b/testing/lin/clagsy.f
@@ -55,7 +55,7 @@
 *
 *  CLAGSY generates a complex symmetric matrix A, by pre- and post-
 *  multiplying a real diagonal matrix D with a random unitary matrix:
-*  A = U*D*U**T. The semi-bandwidth may then be reduced to k by
+*  A = U*D*U^T. The semi-bandwidth may then be reduced to k by
 *  additional unitary transformations.
 *
 *  Arguments
diff --git a/testing/lin/clarhs.f b/testing/lin/clarhs.f
index 4a683f6a72289a29089ea19eef881025ca04fc3e..22165f3a3ca9f9fef93d640ea72ce053c215d131 100644
--- a/testing/lin/clarhs.f
+++ b/testing/lin/clarhs.f
@@ -58,7 +58,7 @@
 *  CLARHS chooses a set of NRHS random solution vectors and sets
 *  up the right hand sides for the linear system
 *     op( A ) * X = B,
-*  where op( A ) may be A, A**T (transpose of A), or A^H (conjugate
+*  where op( A ) may be A, A^T (transpose of A), or A^H (conjugate
 *  transpose of A).
 *
 *  Arguments
@@ -102,7 +102,7 @@
 *          Used only if A is nonsymmetric; specifies the operation
 *          applied to the matrix A.
 *          = 'N':  B := A    * X
-*          = 'T':  B := A**T * X
+*          = 'T':  B := A^T * X
 *          = 'C':  B := A^H * X
 *
 *  M       (input) INTEGER
diff --git a/testing/lin/clatrs.f b/testing/lin/clatrs.f
index 0510829d6f75f906fe424906ca3ffe16771b41ca..38cc5a2a22794de7b5929036b0aaf24be1d165da 100644
--- a/testing/lin/clatrs.f
+++ b/testing/lin/clatrs.f
@@ -57,10 +57,10 @@
 *
 *  CLATRS solves one of the triangular systems
 *
-*     A * x = s*b,  A**T * x = s*b,  or  A^H * x = s*b,
+*     A * x = s*b,  A^T * x = s*b,  or  A^H * x = s*b,
 *
 *  with scaling to prevent overflow.  Here A is an upper or lower
-*  triangular matrix, A**T denotes the transpose of A, A^H denotes the
+*  triangular matrix, A^T denotes the transpose of A, A^H denotes the
 *  conjugate transpose of A, x and b are n-element vectors, and s is a
 *  scaling factor, usually less than or equal to 1, chosen so that the
 *  components of x will be less than the overflow threshold.  If the
@@ -79,7 +79,7 @@
 *  TRANS   (input) CHARACTER*1
 *          Specifies the operation applied to A.
 *          = 'N':  Solve A * x = s*b     (No transpose)
-*          = 'T':  Solve A**T * x = s*b  (Transpose)
+*          = 'T':  Solve A^T * x = s*b  (Transpose)
 *          = 'C':  Solve A^H * x = s*b  (Conjugate transpose)
 *
 *  DIAG    (input) CHARACTER*1
@@ -115,7 +115,7 @@
 *
 *  SCALE   (output) REAL
 *          The scaling factor s for the triangular system
-*             A * x = s*b,  A**T * x = s*b,  or  A^H * x = s*b.
+*             A * x = s*b,  A^T * x = s*b,  or  A^H * x = s*b.
 *          If SCALE = 0, the matrix A is singular or badly scaled, and
 *          the vector x is an exact or approximate solution to A*x = 0.
 *
@@ -181,7 +181,7 @@
 *  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
 *  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
 *
-*  Similarly, a row-wise scheme is used to solve A**T *x = b  or
+*  Similarly, a row-wise scheme is used to solve A^T *x = b  or
 *  A^H *x = b.  The basic algorithm for A upper triangular is
 *
 *       for j = 1, ..., n
@@ -412,7 +412,7 @@
 *
       ELSE
 *
-*        Compute the growth in A**T * x = b  or  A^H * x = b.
+*        Compute the growth in A^T * x = b  or  A^H * x = b.
 *
          IF( UPPER ) THEN
             JFIRST = 1
@@ -632,7 +632,7 @@
 *
          ELSE IF( LSAME( TRANS, 'T' ) ) THEN
 *
-*           Solve A**T * x = b
+*           Solve A^T * x = b
 *
             DO 150 J = JFIRST, JLAST, JINC
 *
@@ -744,7 +744,7 @@
                      ELSE
 *
 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
-*                       scale = 0 and compute a solution to A**T *x = 0.
+*                       scale = 0 and compute a solution to A^T *x = 0.
 *
                         DO 140 I = 1, N
                            X( I ) = ZERO
diff --git a/testing/lin/dpocon.f b/testing/lin/dpocon.f
index 43c957d38dcb564a1de03454ad7b8f0eddf590dc..1a4c1b67aff75cead617c38adba9792e0bd3a4c2 100644
--- a/testing/lin/dpocon.f
+++ b/testing/lin/dpocon.f
@@ -59,7 +59,7 @@
 *
 *  DPOCON estimates the reciprocal of the condition number (in the
 *  1-norm) of a real symmetric positive definite matrix using the
-*  Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF.
+*  Cholesky factorization A = U^T*U or A = L*L^T computed by DPOTRF.
 *
 *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
 *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
@@ -76,7 +76,7 @@
 *
 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
 *          The triangular factor U or L from the Cholesky factorization
-*          A = U**T*U or A = L*L**T, as computed by DPOTRF.
+*          A = U^T*U or A = L*L^T, as computed by DPOTRF.
 *
 *  LDA     (input) INTEGER
 *          The leading dimension of the array A.  LDA >= max(1,N).
diff --git a/testing/lin/dporfs.f b/testing/lin/dporfs.f
index 3a14966380e643851c885f0f8c0223a62e90d0ad..c93d5793ae10f3ff52dca15a531106f49668b68a 100644
--- a/testing/lin/dporfs.f
+++ b/testing/lin/dporfs.f
@@ -92,7 +92,7 @@
 *
 *  AF      (input) DOUBLE PRECISION array, dimension (LDAF,N)
 *          The triangular factor U or L from the Cholesky factorization
-*          A = U**T*U or A = L*L**T, as computed by DPOTRF.
+*          A = U^T*U or A = L*L^T, as computed by DPOTRF.
 *
 *  LDAF    (input) INTEGER
 *          The leading dimension of the array AF.  LDAF >= max(1,N).
diff --git a/testing/lin/dposvx.f b/testing/lin/dposvx.f
index aeca6aee098c20c774b0f430ce3f659d4b85391e..79d723a27095f4205d0fcf03457417ac4f9eaa63 100644
--- a/testing/lin/dposvx.f
+++ b/testing/lin/dposvx.f
@@ -61,7 +61,7 @@
 *  Purpose
 *  =======
 *
-*  DPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
+*  DPOSVX uses the Cholesky factorization A = U^T*U or A = L*L^T to
 *  compute the solution to a real system of linear equations
 *     A * X = B,
 *  where A is an N-by-N symmetric positive definite matrix and X and B
@@ -84,8 +84,8 @@
 *
 *  2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
 *     factor the matrix A (after equilibration if FACT = 'E') as
-*        A = U**T* U,  if UPLO = 'U', or
-*        A = L * L**T,  if UPLO = 'L',
+*        A = U^T* U,  if UPLO = 'U', or
+*        A = L * L^T,  if UPLO = 'L',
 *     where U is an upper triangular matrix and L is a lower triangular
 *     matrix.
 *
@@ -156,18 +156,18 @@
 *  AF      (input or output) DOUBLE PRECISION array, dimension (LDAF,N)
 *          If FACT = 'F', then AF is an input argument and on entry
 *          contains the triangular factor U or L from the Cholesky
-*          factorization A = U**T*U or A = L*L**T, in the same storage
+*          factorization A = U^T*U or A = L*L^T, in the same storage
 *          format as A.  If EQUED .ne. 'N', then AF is the factored form
 *          of the equilibrated matrix diag(S)*A*diag(S).
 *
 *          If FACT = 'N', then AF is an output argument and on exit
 *          returns the triangular factor U or L from the Cholesky
-*          factorization A = U**T*U or A = L*L**T of the original
+*          factorization A = U^T*U or A = L*L^T of the original
 *          matrix A.
 *
 *          If FACT = 'E', then AF is an output argument and on exit
 *          returns the triangular factor U or L from the Cholesky
-*          factorization A = U**T*U or A = L*L**T of the equilibrated
+*          factorization A = U^T*U or A = L*L^T of the equilibrated
 *          matrix A (see the description of A for the form of the
 *          equilibrated matrix).
 *
diff --git a/testing/lin/dpotri.f b/testing/lin/dpotri.f
index f8585b348923865d9cf30ee2b65455ec92c12c55..2a5f4c2dd8fcb52f69cc8a07855cad072e8e2cb4 100644
--- a/testing/lin/dpotri.f
+++ b/testing/lin/dpotri.f
@@ -53,7 +53,7 @@
 *  =======
 *
 *  DPOTRI computes the inverse of a real symmetric positive definite
-*  matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
+*  matrix A using the Cholesky factorization A = U^T*U or A = L*L^T
 *  computed by DPOTRF.
 *
 *  Arguments
@@ -68,7 +68,7 @@
 *
 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 *          On entry, the triangular factor U or L from the Cholesky
-*          factorization A = U**T*U or A = L*L**T, as computed by
+*          factorization A = U^T*U or A = L*L^T, as computed by
 *          DPOTRF.
 *          On exit, the upper or lower triangle of the (symmetric)
 *          inverse of A, overwriting the input factor U or L.
diff --git a/testing/lin/spocon.f b/testing/lin/spocon.f
index 02392607f29753df36ab81a6e9b7180327d0a56c..380896480399b19aeb1ccb731b54e7d5471e11d8 100644
--- a/testing/lin/spocon.f
+++ b/testing/lin/spocon.f
@@ -59,7 +59,7 @@
 *
 *  SPOCON estimates the reciprocal of the condition number (in the 
 *  1-norm) of a real symmetric positive definite matrix using the
-*  Cholesky factorization A = U**T*U or A = L*L**T computed by SPOTRF.
+*  Cholesky factorization A = U^T*U or A = L*L^T computed by SPOTRF.
 *
 *  An estimate is obtained for norm(inv(A)), and the reciprocal of the
 *  condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
@@ -76,7 +76,7 @@
 *
 *  A       (input) REAL array, dimension (LDA,N)
 *          The triangular factor U or L from the Cholesky factorization
-*          A = U**T*U or A = L*L**T, as computed by SPOTRF.
+*          A = U^T*U or A = L*L^T, as computed by SPOTRF.
 *
 *  LDA     (input) INTEGER
 *          The leading dimension of the array A.  LDA >= max(1,N).
diff --git a/testing/lin/sporfs.f b/testing/lin/sporfs.f
index 8dcdea76046fbc0807d4449a3aab69e2ab617e03..e633b09786fb03ae8d62a9dc7f53480700ad89b9 100644
--- a/testing/lin/sporfs.f
+++ b/testing/lin/sporfs.f
@@ -92,7 +92,7 @@
 *
 *  AF      (input) REAL array, dimension (LDAF,N)
 *          The triangular factor U or L from the Cholesky factorization
-*          A = U**T*U or A = L*L**T, as computed by SPOTRF.
+*          A = U^T*U or A = L*L^T, as computed by SPOTRF.
 *
 *  LDAF    (input) INTEGER
 *          The leading dimension of the array AF.  LDAF >= max(1,N).
diff --git a/testing/lin/sposvx.f b/testing/lin/sposvx.f
index b8a94475a4a1e078716c782da412a778087933ee..8a8f53564a41ab269f8016f8e691582e3057bd99 100644
--- a/testing/lin/sposvx.f
+++ b/testing/lin/sposvx.f
@@ -61,7 +61,7 @@
 *  Purpose
 *  =======
 *
-*  SPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
+*  SPOSVX uses the Cholesky factorization A = U^T*U or A = L*L^T to
 *  compute the solution to a real system of linear equations
 *     A * X = B,
 *  where A is an N-by-N symmetric positive definite matrix and X and B
@@ -84,8 +84,8 @@
 *
 *  2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
 *     factor the matrix A (after equilibration if FACT = 'E') as
-*        A = U**T* U,  if UPLO = 'U', or
-*        A = L * L**T,  if UPLO = 'L',
+*        A = U^T* U,  if UPLO = 'U', or
+*        A = L * L^T,  if UPLO = 'L',
 *     where U is an upper triangular matrix and L is a lower triangular
 *     matrix.
 *
@@ -156,18 +156,18 @@
 *  AF      (input or output) REAL array, dimension (LDAF,N)
 *          If FACT = 'F', then AF is an input argument and on entry
 *          contains the triangular factor U or L from the Cholesky
-*          factorization A = U**T*U or A = L*L**T, in the same storage
+*          factorization A = U^T*U or A = L*L^T, in the same storage
 *          format as A.  If EQUED .ne. 'N', then AF is the factored form
 *          of the equilibrated matrix diag(S)*A*diag(S).
 *
 *          If FACT = 'N', then AF is an output argument and on exit
 *          returns the triangular factor U or L from the Cholesky
-*          factorization A = U**T*U or A = L*L**T of the original
+*          factorization A = U^T*U or A = L*L^T of the original
 *          matrix A.
 *
 *          If FACT = 'E', then AF is an output argument and on exit
 *          returns the triangular factor U or L from the Cholesky
-*          factorization A = U**T*U or A = L*L**T of the equilibrated
+*          factorization A = U^T*U or A = L*L^T of the equilibrated
 *          matrix A (see the description of A for the form of the
 *          equilibrated matrix).
 *
diff --git a/testing/lin/spotri.f b/testing/lin/spotri.f
index 13885e2fd868ce9fb6b2d6bd4459839437a62341..d52f05699f40fcb185e124cee9005e390c9db15d 100644
--- a/testing/lin/spotri.f
+++ b/testing/lin/spotri.f
@@ -53,7 +53,7 @@
 *  =======
 *
 *  SPOTRI computes the inverse of a real symmetric positive definite
-*  matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
+*  matrix A using the Cholesky factorization A = U^T*U or A = L*L^T
 *  computed by SPOTRF.
 *
 *  Arguments
@@ -68,7 +68,7 @@
 *
 *  A       (input/output) REAL array, dimension (LDA,N)
 *          On entry, the triangular factor U or L from the Cholesky
-*          factorization A = U**T*U or A = L*L**T, as computed by
+*          factorization A = U^T*U or A = L*L^T, as computed by
 *          SPOTRF.
 *          On exit, the upper or lower triangle of the (symmetric)
 *          inverse of A, overwriting the input factor U or L.
diff --git a/testing/lin/zlagsy.f b/testing/lin/zlagsy.f
index d2a05500dac37b3fccac841b1f89eb5ac300c31b..a9366c90c0bb3eae64e7d50abc109f655d414ad1 100644
--- a/testing/lin/zlagsy.f
+++ b/testing/lin/zlagsy.f
@@ -55,7 +55,7 @@
 *
 *  ZLAGSY generates a complex symmetric matrix A, by pre- and post-
 *  multiplying a real diagonal matrix D with a random unitary matrix:
-*  A = U*D*U**T. The semi-bandwidth may then be reduced to k by
+*  A = U*D*U^T. The semi-bandwidth may then be reduced to k by
 *  additional unitary transformations.
 *
 *  Arguments
diff --git a/testing/lin/zlarhs.f b/testing/lin/zlarhs.f
index 9ace778057818528e04c8d5ff598615bfc4ffba1..333feeb719caa639981cb8dea8404a1208543ce7 100644
--- a/testing/lin/zlarhs.f
+++ b/testing/lin/zlarhs.f
@@ -58,7 +58,7 @@
 *  ZLARHS chooses a set of NRHS random solution vectors and sets
 *  up the right hand sides for the linear system
 *     op( A ) * X = B,
-*  where op( A ) may be A, A**T (transpose of A), or A^H (conjugate
+*  where op( A ) may be A, A^T (transpose of A), or A^H (conjugate
 *  transpose of A).
 *
 *  Arguments
@@ -102,7 +102,7 @@
 *          Used only if A is nonsymmetric; specifies the operation
 *          applied to the matrix A.
 *          = 'N':  B := A    * X
-*          = 'T':  B := A**T * X
+*          = 'T':  B := A^T * X
 *          = 'C':  B := A^H * X
 *
 *  M       (input) INTEGER
diff --git a/testing/lin/zlatrs.f b/testing/lin/zlatrs.f
index 6447bee9aba969bc1f6bebbf2d73ad2c963fa977..ba7f497efb868ac5c9f1411f3228203de223cf56 100644
--- a/testing/lin/zlatrs.f
+++ b/testing/lin/zlatrs.f
@@ -57,10 +57,10 @@
 *
 *  ZLATRS solves one of the triangular systems
 *
-*     A * x = s*b,  A**T * x = s*b,  or  A^H * x = s*b,
+*     A * x = s*b,  A^T * x = s*b,  or  A^H * x = s*b,
 *
 *  with scaling to prevent overflow.  Here A is an upper or lower
-*  triangular matrix, A**T denotes the transpose of A, A^H denotes the
+*  triangular matrix, A^T denotes the transpose of A, A^H denotes the
 *  conjugate transpose of A, x and b are n-element vectors, and s is a
 *  scaling factor, usually less than or equal to 1, chosen so that the
 *  components of x will be less than the overflow threshold.  If the
@@ -79,7 +79,7 @@
 *  TRANS   (input) CHARACTER*1
 *          Specifies the operation applied to A.
 *          = 'N':  Solve A * x = s*b     (No transpose)
-*          = 'T':  Solve A**T * x = s*b  (Transpose)
+*          = 'T':  Solve A^T * x = s*b  (Transpose)
 *          = 'C':  Solve A^H * x = s*b  (Conjugate transpose)
 *
 *  DIAG    (input) CHARACTER*1
@@ -115,7 +115,7 @@
 *
 *  SCALE   (output) DOUBLE PRECISION
 *          The scaling factor s for the triangular system
-*             A * x = s*b,  A**T * x = s*b,  or  A^H * x = s*b.
+*             A * x = s*b,  A^T * x = s*b,  or  A^H * x = s*b.
 *          If SCALE = 0, the matrix A is singular or badly scaled, and
 *          the vector x is an exact or approximate solution to A*x = 0.
 *
@@ -181,7 +181,7 @@
 *  prevent overflow, but if the bound overflows, x is set to 0, x(j) to
 *  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.
 *
-*  Similarly, a row-wise scheme is used to solve A**T *x = b  or
+*  Similarly, a row-wise scheme is used to solve A^T *x = b  or
 *  A^H *x = b.  The basic algorithm for A upper triangular is
 *
 *       for j = 1, ..., n
@@ -412,7 +412,7 @@
 *
       ELSE
 *
-*        Compute the growth in A**T * x = b  or  A^H * x = b.
+*        Compute the growth in A^T * x = b  or  A^H * x = b.
 *
          IF( UPPER ) THEN
             JFIRST = 1
@@ -632,7 +632,7 @@
 *
          ELSE IF( LSAME( TRANS, 'T' ) ) THEN
 *
-*           Solve A**T * x = b
+*           Solve A^T * x = b
 *
             DO 170 J = JFIRST, JLAST, JINC
 *
@@ -744,7 +744,7 @@
                   ELSE
 *
 *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and
-*                       scale = 0 and compute a solution to A**T *x = 0.
+*                       scale = 0 and compute a solution to A^T *x = 0.
 *
                      DO 150 I = 1, N
                         X( I ) = ZERO